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ABSTRACT 

We present a series of 3-D nonlinear simulations of solar-like convection, carried out using the 
Anelastic Spherical Harmonic (ASH) code, that are designed to isolate those processes that drive and 
shape meridional circulations within stellar convection zones. These simulations have been constructed 
so as to span the transition between solar-like differential rotation (fast equator/slow poles) and “anti- 
solar’ differential rotation (slow equator/fast poles). Solar-like states of differential rotation, arising 
when convection is rotationally constrained, are characterized by a very different convective Reynolds 
stress than anti-solar regimes, wherein convection only weakly senses the Coriolis force. We find that 
the angular momentum transport by convective Reynolds stress plays a central role in establishing 
the meridional flow profiles in these simulations. We find that the transition from single-celled to 
multi-celled meridional circulation profiles in strong and weak regimes of rotational constraint is 
linked to a change in the convective Reynolds stress, a clear demonstration of gyroscopic pumping. 
Latitudinal thermal variations differ between these different regimes, with those in the solar-like regime 
conspiring to suppress a single cell of meridional circulation, whereas the cool poles and warm equator 
established in the anti-solar states tend to promote single-celled circulations. Though the convective 
angular momentum transport becomes radially inward at mid-latitudes in anti-solar regimes, it is 
the meridional circulation that is primarily responsible for establishing a rapidly-rotating pole. We 
conclude with a discussion of how these results relate to the Sun, and suggest that the Sun may he 
near the transition between rapidly-rotating and slowly-rotating regimes. 


1. INTRODUCTION 

Mean (longitudinally-averaged) flows play an essential 
role in all recent models of the solar activity cycle. Differ¬ 
ential rotation (the mean zonal flow) generates toroidal 
field from poloidal field and amplifies it, tapping the ki¬ 
netic energy of shear. This is the well-known O-effect and 
it is generally held responsible for the observed dispar¬ 
ity between the toroidal flux inferred from bipolar active 
regions and the much smaller poloidal flux inferred from 
high-latitude magnetograms. 

The role of the mean meridional circulation (MC) in 
the operation of the solar dynamo is more controve rsial 
but potentially no less profound. As discussed in §2.1[ 
mean-field dynamo models have shown that the ampli¬ 
tude and structure of the MC can have a substantial 
influence on cycle properties, and in particular, it may 
largely regulate the cycle period. 

Yet, despite the central role that mean flows play in 
solar and stellar dynamo models, their existence is often 
taken for granted. Most mean-field dynamo models pub¬ 
lished in the literature are kinematic in the sense that 
the velocity field is specified so the origin of mean flows 
lies outside the scope of the model. Non-kinematic mean- 
field models do exist but many take a heuristic approach, 
introducing parameterizations for convective momentum 
and energy transport that produce mean flows w ith pre¬ 
conceived amplitudes and profiles (e.g. Rempel||2006). 

The utility of these kinematic and heuristic approaches 
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to dynamo modeling rests on the availability of suffi¬ 
cient observational data to constrain the mean flows that 
are employed. However, this observational data is lim¬ 
ited. Global helioseismic inversions provide a reliable 
measure of the time-averaged angular velocity profile 
Q,(r,Q) throughout most of the solar convection zone (CZ, 


Thompson et al. 2003 


2.2 1 relatively little is 


Howe 2009), but as discussed in 
mown about the MC. Even less 
is Known about mean flows in other stars (see Q. Yet, 
stars provide an invaluable test of dynamo models; any 
viable model of the solar activity cycle must also account 
for cyclic magnetic activity in solar-like stars. 

Without some theoretical guidance as to what merid¬ 
ional and zonal flows one should expect to find in solar 
and stellar convection zones, kinematic dynamo models 
cannot completely account for observed patterns in solar 
and stellar magnetic activity. Furthermore, theoretical 
guidance is needed to help motivate and interpret ongo¬ 
ing observat iona l investigations of the solar meridional 
circulation (£2.2) which seek to probe more generally the 
internal dynamics of the Sun and the origins of solar 

magnetic activity. __ 

In a previous paper, (Miesch et al.||2012| hereafter Pa¬ 
per 1), we sought such theoretical guidance based only 
on the observed structure of solar mean flows and ro¬ 
bust dynamical balances obtained from the fundamen¬ 
tal magnetohydro dynamic (MHD) equations. From this 
alone we were able to identify gyroscopic pumping as 
the principle physical mechanism that maintains the so¬ 
lar MC. This work is summarized in [J2) However, fur- 
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ther progress requires an understanding of turbulent so¬ 
lar convection and in particular the convective angular 
momentum transport. To this end, in this paper, we 
have initiated a series of global convection simulations 
designed to elucidate the nature of convective transport 
in rotating, stratified, spherical shells and the mean flows 
it produces. We focus here on solar-like stars with con¬ 
vective envelopes but similar dynamics are likely to occur 
in other types of stars with different convection zone ge¬ 
ometries. 

In ^2]we survey how this series of simulations fits within 
the context of previous theoretical and observational 
work and in ^3] we discuss the motivation and method¬ 
ology behind how the numerical experiments were set 
up. We then present an overview of our results in Q 
identifying two distinct dynamical regimes that exhibit 
qualitati vely different mean flows. As in previous stud¬ 
ies (see j |2.3| ), the dynamical regimes are distinguished 
by the degree of rotational influence on the convection, 
as quantified by the Rossby number Ro = U/(2CIqL), 
where U and L are characteristic convective velocity and 
length scales and is mean rotation rate. Fast rotators 
produce solar-like H profiles (fast equator, slow poles) 
with multi-celled MC profiles while slow rotators pro¬ 
duce anti-solar D profiles (slow equator, fast poles) with 
single-celled MC profiles (poleward flow in the upper CZ 
and equatorward in the lower CZ). We demonstrate that 
these slow and fast regimes can be achieved not only by 
changing the rotation rate but also by changing the me¬ 
chanical and thermal dissipation, which in turn impacts 
the buoyancy driving. In ^5] we interpret the resulting 
MC profiles as a consequence of the changing nature of 
the convective Reynolds stress, within the context of gy¬ 
roscopic pumping. We address the establishment and 
significance of thermal gradients in © 

Which regime is the Sun in? Since it has a solar-like 
O profile by definition, one would expect it to be classi¬ 
fied as a rapid rotator. However, the situation is more 
subtle than this. In © we demonstrate that the tran¬ 
sition between solar and anti-solar differential rotation 
is intimately linked to the MC and may be sensitive to 
the complex dynamics occurring in the boundary layers, 
namely the tachocline and the near-surface shear layer. 
The Sun lies very close to this transition, which may 
imply that it has somewhat atypical mean flow profiles. 
This in turn has important implications for solar dynamo 
models and their application to other stars. These issues 
too are addressed in Sj7]and our principle results and con¬ 
clusions are summarized in ^8j 


2. THEORETICAL AND OBSERVATIONAL BACKGROUND 
2.1. Meridional Circulation and the Solar Dynamo 

The potential significance of the MC to the opera¬ 
tion of the solar dynamo is apparent from the observed 
evolution of the radial magnetic flux in the solar pho¬ 
tosphere. In particular, the cyclic polarity reversals of 
polar fields observed at the solar surface are often at¬ 
tributed to the poleward transport of the trailing mag¬ 
netic flux in emerging active regions due to the com¬ 
bined action of the meridional flow and turbulent dif¬ 
fusion dSh eeley 2005 Baumann et al.||2006; Schrijver & 
Liu 2008|). Though the timing of the most recen t north 
polar reversal in 2012 challenges this paradigm (Shiota 


et al. 2012), it is clear that the meridional flow plays an 
important role in redistributing the magnetic flux that 
threads through the solar surface. 

These observations and models of surface flux trans¬ 
port help to justify the prominent role that the MC plays 
in many recent dynamo models of the solar cycle. This 
role is most dramatic in a class of models known as flux- 
transport (FT) dynamo models, in which the equator- 
ward advection of toroidal flux by the meridional flow 
in the CZ and tachocline accounts for the observed mi¬ 
gration of sunspots toward the equator over the co urse 
of a cycle, know n as the solar butterfly diagram 
fc Sheeley] |1991| |Dikpati fc Charbonneau|| 19991 P^ uker 
et al. 20011 Dikpati fe Gilman||20061 |Yeates et al. 


Munoz-Jaramillo et al.||2009[ Jouve & Brun|2007|). By 
regulating the solar butterfly diagram, the MU also reg¬ 
ulates the cycle period. 

In a subset of flux-transport dynamo models that op¬ 
erate in the so-called advection-dominated regime, the 
MC plays another important role. In practice, most ex¬ 
isting advection-dominated flux-transport (FT) dynamo 
models are also Babcock-Leighton (BL) dynamo models 
in which the main source of poloidal flux lies in the buoy¬ 
ant destabilizatio n, emergence and subsequent dispersal 
of tor oidal flux (|Dikpati & Gilman 2009| Charbonneau 
2010). Mean-field parameterizations of the BL mecha- 


nism typically take the form of a polo i dal so ur ce te rm 
confined to the surface layers (iDikpati fe Charbonneau 


1999 Rempel 2006 Munoz-Jaramillo et al. 2010[ e.g.J. 
Thus, in order Tor the dynamo to operate, there must 
be a coupling mechanism that transports magnetic flux 
from the poloidal source region near the surface to the 
region where the toroidal field is generated, which gener¬ 
ally lies in the lower CZ. In advection-dominated BL/FT 
dynamo models, this coupling mechanism is the MC. 

Changes in the amplitude and structure of the merid¬ 
ional flow in FT dynamo models can have a dramatic 
effect on the length and amplitude of magnetic cycles. 
Faster meridional flows or multiple cells in latitude can 
shorten the cycle period and influence the polar field 
strength while slower flows can lengthen cycles and sup¬ 
press activity, potentially induc i ng a grand minimum 


press activity, potentially induc i ng a grand mn 
(jDikpati fc Charbonneau 19991 Jouve fc Brun 


bchrijver & Liu'2008 Jiang et al. 


Karak 


2007 


2010 ) 


Multiple cells in radius can have an even greater impact, 
disrupting the operation of the dynamo and p otentially 
rever sing the sense of the butterfly diagram (Jouve & 
Brun 2007). 

The profound role of the MC in contemporary solar 
dynamo models has motivated intense, ongoing observa¬ 
tional efforts to probe the meridional flow structure in the 
deep CZ. However, this is a formidable challenge, largely 
because the flow is expected to be so weak; typical am¬ 
plitudes used in FT dynamo models are only 2-3 ms -1 
near the base of the CZ. Thus, observational constraints 
are currently limited but may improve with more data, 
as we now address. 

2.2. Observational Context 

One of the most notable triumphs of helioseismology 
has been the mapping of the internal solar rotation profile 
9 as a function of radiu s r and colatitude 6 (Thompson 
et al.||2003 Howe||2009). These helioseismic rotation in¬ 
versions span most of the solar convection zone, though 
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they become less reliable near the poles and in the deep 
interior due both to the coarser resolution of the inver¬ 
sion kernels and the smaller amplitude of the intrinsic 
frequency splitting near the rotation axis. They reveal a 
monotonic decrease of D by about 30% from equator to 
pole that persists throughout most of the CZ, with lit¬ 
tle radial variation apart from two boundary layers near 
the bottom and top, known as the tachocline and the 
near-surface shear layer respectively. 

The most well-established result for the solar MC, 
meanwhile, is the existence of a systematic poleward flow 
of 10-20 ms -1 at the surface at latitudes below about db 
50-60°. However, different techniques often give different 
results as to the subsurface structure, the high-latitude 
structure, and even the low-latitude amplitude. For ex¬ 
ample, measurements of the meridional flow based on 
the correlation tracking of magnetic features generally 
give somewhat lower flow speeds than those ba sed o n 
Dopp l er shifts and local hel ioseismic inversions (Ulrich 
2010). Dikpati et al. (2010b) attribute this difference to 
the effects of turbulent diffusion. Much of the variation in 
amplitude estimates arises from the intrinsic variability 
of the flow itself, which is known to change with the solar 


cycle ([Ul rich 2010" Basu fc Antia||2010 Hathaway||2010 
Hathaway & Rightmire||201l|). However, these solar cy¬ 
cle variations can be understood as perturbations about 
a persistent background circulation (poleward near the 
surface at low-mid latitudes) that are induced by surface 
magnetism and its impact on radiative cooling. In this 
paper we are interested only in the persistent background 
circulation so we disregard the solar cycle variations. 

At higher latitudes t he res ul ts ar e less clear. The recent 
Doppler analysis by Ulrich (2010) suggests the possible 
presence of a high-latitudecounter-cell with equatorward 
surface flow at latitudes above about 60°. This feature 
is most apparent ne ar sol ar minimum and may be in¬ 
termittent. Indeed, Dikpati et al. (2010a) argue that 
the absence of such a high-latitude counter cell through 
much of solar cycle 23 may account for the relatively 
deep minimum of 2008-2009 and the delayed onset of so¬ 
lar cycle 24. However, evidence for such a high-latitude 
counter-cell based on feature tracking and helioseismic 
inversions has so far been lacking, with the former sug¬ 
gesting that pol eward flow may pers ist all the way to the 
poles (Rightmire-Upton et al.l 2012t and the latter being 
thus far inconclusive (|Zhao et al.||20f2j). Further insight 
into this question should be possible in the coming years 
with the availability of longer time series from the Helio¬ 
seismic Magnetic Imager (HMI) on the Solar Dynamics 
Observatory (SDO), which has the higher spatial resolu¬ 
tion needed to minimize limb effects, and also with the 
advent of future missions such as Solar Orbiter, which 
will have a high-latitude vantage point. 

Of these three principle techniques for measuring the 
meridional flow (Doppler measurements, feature track¬ 
ing, and local helioseismic inversions), only helioseismic 
inversions and feature tracking can probe the meridional 
flow below the surface. While the observational consen¬ 
sus is that there is a systematic poleward flow of about 
10-20 ms -1 at latitudes below about 60° in the near¬ 
surface layers of the Sun, the subsurface structure is more 
controversial. The fir st s tudies based on SOHO/MDI 


ies suggested that this poleward flow may extend signif¬ 
icant lyMnepe^smmm^ a^JeasWthe upper half of the 
CZ (Braun & Fan 1998 Chou & Dai 2001). However, 
more recent observations and inversions based on longer 
time series, new analysis techniques, and new data from 
H MI/SDO have begun to cast doubt on this conclusion. 


Hathaway (2012b) argued based on surface feature 
tracking that the return equatorward flow required for 
mass conservation may occur as shallow as 0.93R. In¬ 
ferring subsurface flows with feature tracking relies on 
the hypothesis that larger horizontal scales in the pho¬ 
tosphere trace convective structures that extend deeper 
into t he CZ and may th u s be advected by de eper mean 
flows ( Hathaway||2012a ). Hathaway (2012a|b) tests this 
hypothesis by comparing interred zonal flows to helio¬ 
seismic rotational inversions and then applies it to infer 
the meridional flow. Results suggest a flow reversal from 
poleward to equatorward at a depth of roughly 50 Mm 
(r ^ 0.93, where R is the solar radius). This is notably 
shallow compared to earlier estimates obtained from lo¬ 
cal helioseismic inversions. 

More recent helioseismic investigations, fueled by the 
high-quality data from HMI/SDO and by the availability 
of long time series from GONG and SOHO/MD^ sug¬ 
gest that this picture ma y be changi ng. Time-distance 
helioseismic inversions by Zhao et al. ( 2013]) , taking into 
account systematic center-to-limb effects on the inver¬ 
sions, appear to support this conclusion, showing evi¬ 
dence for multiple reversals in the latitudinal flow com¬ 
pon ent, th e shallowest of w hich being 0.91 R (see also 
Mitra-Kraev & Thompson 2007|. However, such a shal¬ 
low reversal is not found in global inversions of low-order 
oscillation modes w here the meridional flow enters as 
higher-order effect (|Schad et al.||2012|). 

Given the high stakes from the perspective of solar 
and stellar dynamo theory, the current uncertainty with 
regard to the meridional flow structure in the deep con¬ 
vection zone, and the vibrancy and promise of ongoing 
observational efforts, further theoretical insight is sorely 
needed. Here we seek such insight from convection sim¬ 
ulations, after a brief overview of the fundamental pro- 
involved (©• 


cesses 


2.3. Dynamical Models 

As noted above, kinematic mean-field dynamo models 
cannot provide any insight into the origin and structure 
of the MC. For this, dynamical models are required. Such 
models must address the transport of energy and momen¬ 
tum by convection and the subsequent nonlinear feed¬ 
backs involving th e MC, differential rotati on, and ther¬ 
mal stratification (Miesch & Toomre 2009). This can be 
done either explicitly through global 3-D convection sim¬ 
ulations or implicitly through non-kinematic mean-field 
models that include parameterizations for the turbulent 
transport. 

Global convection simulations have revealed two dis¬ 
tinct dynamical regimes characterized by q ualitatively 


19771 Glatzmaier & Gilman 1982| 

jAurnou et al.[ 

2007 

Brun & Palacios 2009 Matt 

et al.| 

201T Gastine et al. 


data by Giles et al. (1997) suggested that the poleward 
flow extended down to at least 0.9 6R. Subsequent stud- 


1 the Global Oscillations Network Group (GONG) and the 
Michelson Doppler Imager (MDI) on the Solar and Heliospheric 
Observatory (SOHO) 
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20131 Guerrero et al.||2013a| |Gastine et aL|[20l4| |Kapyla| 


et al. 2014 K At strong rotational influence Ro < 1 

the differential rotation is “solar-like” in the sense that 
the equator spins faster than the poles. In other words 
AD = Vt eq — Q p > 0 where Q eq is the equatorial rotation 
rate at the outer boundary of the domain, and Q p is some 
measure of the rotation rate at high latitudes, poleward 
of 70°. For the weakly rotating regime, Ro > 1, the dif¬ 
ferential rotation is referred to as “anti-solar”, meaning 
AD < 0. The anti-solar regime has often been attributed 
to the tendency for the convection to conserve angular 
momentum lo cally, spinning up as motions appr oach the 
rotation axis (|Gilman||1977| |Aurnou et al.||2007|) . 

Though these regimes can be achieved by varying 
the rotation rate D (or its non-dimensional equivalent), 
several studies have demonstrated that the transition 
from solar to anti-solar differential rotation can also be 
achieved by keeping D (or its non-dimensional equiva¬ 
lent) fixed and changi ng the convective driving through 
the R ayleigh number (|Gastine et al.||2014 Kapyla et al 
|2014|). This tends to increase Ro by decreasing the con- 
vective turnover time L/U. In other words, convective 
velocity scales increase while lengt h scales decrease. Near 
the solar/ant i -solar transition, |Gastine et al.| (|2014|) and 
Kapyla et al. (2014) also find that the DR exhibits bista¬ 
bility, with Doth solar and anti-solar regimes possible de¬ 
pending on initial conditions and hysteresis. 

Many of these studies reported a qualitativ e ch ange in 
the MC profile in the two regimes (iMatt et al. 12011; Guer- 

Kap^l 


Gastine et al .|2 013, STl4; iE Lapyla et al 


rero et al. 


2014p. The anti-solar regime (Ro < 1) often exhibits 

a single circulation cell per hemisphere, with poleward 
flow in the upper CZ and equatorward flow in the lower 
CZ. Meanwhile, the solar regime (Ro > 1) often exhibits 
multiple cells in latitude and radius, aligned with the 
rotation axis outside the tangent cylindei^ 

Recent high-resolution models have alsobegun to cap¬ 
ture aspects of the solar near-surface shear layer (NSSL) 
wherein small-scale convection i n the surface layers es¬ 
tablishes an inwa rd 11 gradient (Guerrero et al. 2 013b[ 
Hotta et al. 2014). As suggested oy~ Miesch & Hindman 
(ponjr the JNSSL in these simulations is established by 
a radially inward angular momentum transport, is main¬ 
tained by meridional Reynolds stress, and is intimately 
linked to a polew ard meridional flow. This will be dis¬ 
cussed further in §7.2| below. 

The solar and anti-solar regimes in global convection 
simulations stand in stark contrast to mean-field models 


which typically produce solar-l i ke DR profiles (e.g. Kiiker 
et al. 201 1| Kitchatinov||2012| |Kitchatinov fc Olemskoy 
2012|). Anti-solar profiles generally only occur if there 
Ts a strong, baroclinic, single-celled meridional flow in¬ 
duced by additional factors such as pol ar starspots or 
tidal forcing from a binary companion (|K itchatinov fc 
Riidiger||2004| |Riidiger & Kitchatinov||2007|) . 

Mean-held models of solar-like stars generally exhibit 
single-celled MC profiles but their origin and structure 
vary depending on the modeling approach. One ap¬ 
proach is to seek steady-state solutions of the mean- 
field equations with heuristic or model-based prescrip¬ 
tions for the convective Reynolds stress and heat trans¬ 


port. In such models the MC profile is largely deter¬ 
mined by departures from thermal wind balance due to 


eled as a turbulent diffusion (Kitchatinov 

2012 

Dikpati 

2014). This can lead to nonlocal driving oi 

: the MG from 


or to hig h-latitude cou nter-cells for steep density stratiti- 
cations (Dikpati 2014). However, in this case, the form of 
the MC profile is sensitive to the nature of the turbulent 
viscosity profile that is imposed. This alone breaks the 
degeneracy of the thermal wind b alanc e equation with 
respect to the meridional flow (see j |5.1| below). 

When the full time-dependent mean-held equations are 
solved, the driving of the meridional flow is fundamen¬ 
tally different. Here the MC profile is determined mainly 
by the non-diffusive component of the convective angu¬ 
lar momentum transport, which is typically mo deled b y 


the A -effect in the zonal momentum equation (Rempel 


2005 2007). This is the process of gyro scopic pumping 
that will be described in detail in £5.1 and that is re¬ 
sponsible for the establishment of MC in our convection 


simulations, as demonstrated in £5.2 


3. DESCRIPTION OF NUMERICAL CONVECTION 
EXPERIMENTS 

We present a series of numerical experiments designed 
to examine solar-like convection and its resulting mean 
flows under conditions that span a range of rotational 
constraints and convective efficiency. We model the so¬ 
lar convection zone using the anelastic approximation 
(Gough 1969; Gilman & Glatzmaier 1980). Such an 
approach is appropriate in deep stellar interiors where 
plasma motions are subsonic and perturbations to ther¬ 
modynamic variables are small compared to their mean, 
horizontally averaged values. As thermodynamic per¬ 
turbations are small, thermodynamic variables are lin¬ 
earized about a spherically symmetric reference state 
with density p, pressure P, temperature T. and specific 
entropy S. Fluctuations about this state are denoted as 
p, P, T, and S. In the uniformly-rotating reference frame 
of the star, the anelastic equations under these assump¬ 
tions may be expressed as 


V • (pv) = 0, (1) 


P 


Dv 


+ 2 D 0 x v 


and 


-\7P + pg-\7-V, (2) 


- DS 

pT— =V • [Kprvs + ft r pc p V(T + T)} 


2 pu 




- o (V-U 


(3) 


The velocity v expressed in spherical coordinates is 
v = (u r ,U 0 ,u^) relative to a frame rotating at constant 
angular velocity D 0 , g is the gravitational acceleration, 
and V is the viscous stress tensor given by 


2 The tangent cylinder is a cylindrical surface parallel to the 
rotation axis and tangent to the base of the CZ 


T>ij = -2 pv 


e ij ~ ~(V • v)6i 


(4) 
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where e^- is the strain rate tensor. The kinematic viscos¬ 
ity is denoted by v and the thermal diffusivity by k. The 
radiative diffusivity is indicated by ft r , c p is the specific 
heat at constant pressure, and D/Dt is the advective 
derivative. K r and c p remain constant in time. 

This set of equations is closed by assuming the ther¬ 
modynamic fluctuations satisfy the linear relations 


p_P_T_F^_S^ 
p P T 7 P c p 

assuming the ideal gas law 

P = npf , 


(5) 

(6) 


where 7Z is the gas constant. 

We evolve the anelastic equations using the anelastic 
spherical harmonic (ASH) code, details of which may be 
found in Clune et al. (1999) and Brun et al. (2004). ASH 
solves the 3-D MHD equations in a rotating, spherical 
shell using a pseudo-spectral approach. Spherical har¬ 
monics are employed in the horizontal dimension, and 
4th-order finite-differences in the radial direction. ASH 
employs a poloidal/toroidal decomposition for the mass 
flux vector to satisfy the divergence-free constraint of eq. 
[TJ with 


pv = V x V x (We r ) + V x (Ze r ). 


(7) 


Here W and Z are the poloidal and toroidal stream- 
functions respectively, and are the quantities that are 
explicitly evolved in time (as opposed to the individual 
components of v). The unit vector in the radial direction 
is indicated by e r . 

The suite of convection zone models that we exam¬ 
ine are identical in all respects save for the diffusivities 
and the rotation rate. The reference state temperature, 
pressure, and functional form of K r are d erived fro m a 
one-dimensional solar structure model (Brun et al. 2002). 
While the specific heat c p of the one-dimensional struc- 
ture model varies weakly throughout the solar convec¬ 
tion zone, we have adopted a constant value of c p , taken 
from the base of that model. This requires us to itera¬ 
tively solve for the reference state density p, using the 1 - 
D model density as an initial guess, until a density profile 
that satisfies both hydrostatic balance and our equation 
of state is reached. Our simulated convection zones span 
from 0.72 R 0 to 0.965 R 0 , extending to 25 Mm below 
the solar surface, and encompassing roughly 3.5 density 
scale heights. Convective perturbations about this refer¬ 
ence state remain small throughout the simulation, and 
thus we choose not to update the reference state, holding 
it fixed throughout each calculation. 

As the microscopic viscous and thermal diffusivities 
appropriate for plasma motions within the convection 
zone are beyond the capabilities of current computational 
models, our models should be viewed as large eddy sim¬ 
ulations. We adopt a radial profile for v and k with each 
quantity diminishing in depth as v ~ p~ 1 ' 2 . Simulations 
are thus more diffusive near the surface where we ex¬ 
pect the unresolved scales of convective motion to play 
a prominent role in the transport of heat and momen¬ 
tum. This formulation does not, however, capture any 
non-dissipative behavior of the unresolved scales. 

Impenetrable and stress free boundary conditions are 


adopted for each simulation such that 


d{v e /r) _ div^/r) 


dr 


dr 


0|r« 


r bot1Ttop ' 


(8) 


The radial entropy gradient is forced to vanish at the 
lower boundary of the convection zone, and the entropy 
perturbations are forced to vanish at the upper boundary, 
with 

dS 

- 7 ^ = 0 | r=r bot , S = 0\ r=rtop . (9) 


Thus, there is no diffusive entropy flux across the lower 
boundary. Heat is transported across this boundary by 
the radiative heat flux instead, oc ft r VT, which drops to 
near zero at the upper boundary (see Fig. [To| below). The 
vanishing of the entropy gradient at the lower boundary 
implies that this is, by definition, the base of the CZ, 
with no overshoot region. The entropy gradient steepens 
toward the top of the domain, where most of the convec¬ 
tive driving occurs. This will be discussed further in ^ 6 j 
With the value of entropy pinned at the upper boundary, 
the entropy gradient there is free to vary with time. As 
the simulation equilibrates, this gradient attains a sta¬ 
tistically steady value that is sufficient to expel one solar 
luminosity through the upper boundary via thermal dif¬ 
fusion. 

A summary of our input parameter space for those 
cases rotating at the solar rate of S4 0 =2.6xlO -6 s -1 is 
provided in Table [3j We vary the viscosities by a fac¬ 
tor of four and the thermal diffusivities by a factor of 
two for these simulations. The naming convention is 
such that the letter denotes the level of thermal diffu¬ 
sion, with A corresponding to the highest thermal diffu¬ 
sivity of 3.2xlO 13 cm 2 s -1 and C corresponding to the 
lowest level (1.6xlO 13 cm s _1 ). The numbers 0-3 indi¬ 
cate the level of viscous diffusion, with 0 corresponding 
to the lowest kinematic viscosity ( 2 xl 0 12 cm s -1 ), and 
3 corresponding to the highest value of 8xl0 12 cm s _1 . 
Case B 2 , which lies in the center of our parameter space, 
serves as the point around which rotation rate is varied 
for the results presented in 0 We have also provided 
a summary of the resulting energy balances, as well as 
their relative ratios, for these cases in Table 2. There, 
we report on the globally-averaged (over the full spher¬ 
ical shell) convective kinetic energy density (CKE), the 
energy associated with the differential rotation (DRKE), 
and the energy associated with the meridional circulation 
(MCKE). We define these as 


CKE = h [(t; r - (tv )) 2 + iy e - (tv )) 2 + (v# - (tv)) 2 ] , 

( 10 ) 

MCKE = l -~p [(tv ) 2 + (tv) 2 ] , (11) 

and 

DRKE = l -p [( v 0 ) 2 ] , (12) 


where the brackets denote averages over longitude. 

Each model discussed in this paper has been evolved 
for several thousand days, with elapsed simulated time 
ranging from 7,000 days to 14,000 days depending on 
the individual case. This corresponds to several hundred 
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TABLE 1 
Variable Input Parameters 


Case 

Vtop 

fttop 

Pr 

Ro 

Ra/10 5 

IV 

Imax 

Al 

4 

32 

1/8 

0.15 

1.79 

2.61 

340 

A2 

6 

32 

3/16 

0.13 

1.19 

2.61 

170 

A3 

8 

32 

1/4 

0.11 

0.89 

2.61 

85 

BO 

2 

24 

1/12 

0.16 

6.35 

3.48 

680 

BO.5 

3 

24 

1/8 

0.16 

4.24 

3.48 

340 

B1 

4 

24 

1/6 

0.15 

3.18 

3.48 

340 

B2 

6 

24 

1/4 

0.13 

2.12 

3.48 

170 

B3 

8 

24 

1/3 

0.12 

1.59 

3.48 

170 

Cl 

4 

16 

1/4 

0.24 

7.14 

5.22 

340 

C2 

6 

16 

3/8 

0.20 

4.76 

5.22 

340 

C3 

8 

16 

1/2 

0.12 

3.57 

5.22 

170 


Note. — Input parameters for all cases. Viscous and 
thermal diffusivities (y and ft respectively) are quoted in 

n viitn /^v-P 1 ^ n ^ f I 1 L r\ L) -v» n d fl -nil w>V\n-p Dv* in 


4.1. Identification of Mean Flow Regimes 

We find that the structure of meridional circulation re¬ 
alized in our convection zone models is intimately linked 
to the differential rotation established. This relationship 
is most readily illustrated through a series of experiments 
wherein the stellar rotation rate is varied. Using Case 
B2 as our starting point, we have explored rotation rates 
ranging from 2U© to 0.75£2©. 

We plot the resulting differential rotation (upper row ) 
and meridional circulation (lower row ) profiles for this se¬ 
ries of simulations in Fig. [l] At the higher rotation rates, 
most notably 2U© and U©, the resulting differential ro¬ 
tation is decidedly solar-li ke, p ossessing a fast equator 
and retrograde poles (see £2.3). We note that contours 


in these simulations are noticeably cylindrical with re¬ 
spect to the Sun. This may be attributed to the lack of 
an underlying stable zone which is thought to promote 


stant with depth. The Rossby number Ro= v/(2QL) is 

pel 2005; 

Balbus & Schaan 

2012 

). At rotation rates less 

quoted at mid-convection zone. The Rayleigh number Ra= 

than U©, 

the solar-like dif 

ferential rotation transitions 


(dp/dS)(dS/dr)gL 4 /(pnn) is quoted at the top of the sim¬ 
ulation. The length scale L is taken to be the depth of the 
convection zone (1.72xl0 10 cm) in both cases. The modi¬ 


fied Rayleigh Number IV — g \dS c /dr\ /(c p Q 2 ) 
mid-convection zone. 


TABLE 2 
Energy Balance 


is quoted at 


Case 

CKE 

DRKE 

DRKE/CKE 

MCKE 

MCKE/CKE 

Al 

2.09 

2.94 

1.4091 

0.019 

0.0091 

A2 

1.46 

2.38 

1.6241 

0.012 

0.0079 

A3 

1.22 

1.77 

1.4447 

0.010 

0.0080 

B0 

4.39 

0.59 

0.1345 

0.032 

0.0073 

B0.5 

3.52 

1.87 

0.5305 

0.026 

0.0075 

B1 

2.99 

1.73 

0.5777 

0.024 

0.0079 

B2 

2.23 

1.79 

0.8024 

0.017 

0.0078 

B3 

1.96 

1.31 

0.6721 

0.016 

0.0084 

Cl 

4.66 

9.92 

2.1285 

0.054 

0.0115 

C2 

3.73 

6.04 

1.6202 

0.044 

0.0117 

C3 

2.90 

0.18 

0.0625 

0.025 

0.0085 


Note. — Globally averaged energy balance as realized in 
each run. Energy densities have been averaged over 10 rota¬ 
tion periods at the end of each simulation, and are reported 
in units of 10 6 erg cm -3 

rotational and convective overturning times. All simu¬ 
lations have been well equilibrated viscously and ther¬ 
mally. As Prandtl numbers are always less than unity, 
the viscous diffusion timescale is the most limiting of 
these two timescales. With the exception of case BO, 
which has been run for one viscous diffusion time, all 
simulations have been run for at least two viscous dif¬ 
fusion times. We find that case BO is transitioning to 
an anti-solar differential rotation rate after one viscous 
diffusion time of evolution. As we are able to explore 
anti-solar regimes adequately with the more computa¬ 
tionally tractable cases Cl and C2, we have chosen not 
to pursue the evolution of BO further, but include it in 
Tables 1 and 2, as well as Figure [3j to help delineate the 
point of transition between solar and anti-solar behavior. 


4. 


DISTINCT MEAN FLOW REGIMES FOR FAST AND 
SLOW ROTATORS 


to a state characterized by a weakly retrograde equator 
and rapidly prograde poles. These prograde polar states 
tend to possess a much stronger equator-to-pole rota¬ 
tional contrast than their solar-like counterparts. Simi¬ 
lar results were also reported b y other authors , includ¬ 
ing 


Gastine et al. 


Kapyla et al. ( |2U14fr 


.! (!2013| [2014|, 
l# - 


Guerrero et al. (2013a) 


As the lower"row of Fig. [l] is traversed from left to 
right, circulatory patterns assume two distinct morpholo¬ 
gies and exhibit a clear correlation with the structure of 
the differential rotation. In the fast-equator cases, the 
low latitudes are characterized by a counter cell of cir¬ 
culation occurring near mid-convection zone. High lati¬ 
tudes in these cases are largely characterized by a single 
cell of circulation (poleward at the surface), with some 
hints of multi-cellular behavior at high latitudes for the 
case with the fastest rotation. As the rotation rate of the 
system decreases, the high latitude cell extends to lower 
latitudes, ultimately supplanting the counter cell in the 
equatorial regions for rotation rates lower than U©. The 
shift from multi-cellular MC profiles in the solar regime 
to single-celled profiles in th e anti-solar regime has also 
been noted by other authors (|Guerr ero et al |2013 
tine et ak||2014 Kapyla et al.||2(J14p . 


Gas- 


4.2. Transition Between Mean Flow Reqimes at Fixed 

The series of simulations illustrated in Fig. [l] possess 
a common mid-core rms velocity amplitude of about 100 
ms -1 . The transition in mean-flow regimes spanned by 
those simulations is thus most naturally characterized 
by a transition in timescale ratios, namely the ratio of 
the convective overturning time to the rotation period. 
This is quantified by the Rossby number, introduced in 
|l] We can similarly span this transition by holding ft 
fixed and varying the amplitude of the convective veloc¬ 
ities. For systems with a fixed luminosity, such as those 
under consideration here, this is most naturally accom¬ 
plished by varying the transport coefficients v and ft, and 
through them the Rayleigh number Ra. In so doing, we 
can also explore how our transport parameters v and 
ft contribute to the balances realized in the slowly- and 
rapidly-rotating regimes by sampling a range of resulting 
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Fig. 1. — (a-e) Differential rotation profiles (Q — 12 frame ) f° r cases rotating at different fractions of the solar rotation rate 12o (indicated) 
achieved by spinning up and spinning down case B2, which rotates at the solar rate. Rotation rate decreases from left to right, and units are 
indicated in nHz. (f-j) Meridional circulation profiles corresponding to each of the differential rotation profiles in the upper row. Red tones 
denote counter-clockwise motion, and blue tones, clockwise. Profiles have been averaged in time over roughly 200 days in each instance. As 
the Sun is spun down, the polar regions develop prograde rotation, and the equatorial regions retrograde. Similarly, multi-cellular profiles 
of meridional circulation transition into single-celled profiles at low rotation rates. 


Reynolds numbers. 

Others have shown that the transition between the so¬ 
lar and anti-solar rotation regimes can be achieved by fix¬ 
ing the rotation rate and changing the convective driving 
by varying the Rayleigh number (Gilman 1977| Aurnou 

contrast in models such as ours, analogous to the temper¬ 
ature contrast in RBC systems, is not held fixed, but is 
an output quantity, much as the Reynolds number. This 
entropy contrast is built entirely in the upper boundary 
layer, resulting from the entropy gradient realized there 

et al. 

2007; Guerrero et al.||2013a| Gastine et aT 2011; 

and the width of that boundary layer. Convective driv¬ 

Kapy 

a et al. 2014). The key parameter in determin- 

ing is thus accomplished primarily in the upper boundary 

ing the transition is the Rossby number. However, since 

layers of these systems, a situation more analogous to the 


this is an output of the simulation rather than an input, 
the transition is often described in terms of the modified 
Rayleigh number 77* (see Table [3]), which is a measure 
of the re lative streng ths the buoyancy and Coriolis force 
(Gastine et al. 2013). The simulations presented in this 
section will be used m sections [5] to explore the structure 
and origin of the meridional circulation, a topic that has 
received less attention in previous studies. While our 
demonstration of the solar/anti-solar transition is itself 
not new, we find it useful to revisit some aspects of this 
transition as they also bear on the resulting meridional 
circulation. 

It is worth mentioning that the Rayleigh number in 
systems such as these, where entropy is held fixed at one 
boundary, and the entropy gradient fixed at the other, 
possesses an implicit dependence on k not found in clas¬ 
sical Rayleigh-Benard convection (RBC). The entropy 


Sun than RBC flow where driving occurs at the upper 
and lower boundaries. 

While entropy is fixed at the top, the entropy gradient 
is implicitly determined by the solar luminosity, being 
forced to satisfy, in a time-averaged sense, the relation 

^1'— - dt <i3 » 

We thus choose to define the our Rayleigh numbers using 
ds/dr at the upper boundary because this quantity is a 
well-defined simulation input. The consequence is that 
our Rayleigh number now varies as instead of the 
more traditional kT 1 . 

The variation of differential rotation for each case is 
illustrated in Fig. [2| There, profiles of angular velocity, 
averaged over several overturning times, are laid out in 
a grid, with k increasing along the vertical axis and v 
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Fig. 2. — Differential rotation profiles (Q — Q©) for the a,b , and c series of simulations. Case names are indicated and are arranged such 
that viscosity v increases from left to right, and thermal diffusivity k, increases from bottom to top. The most laminar case (A3) appears 
in the upper right, and the most turbulent case (Cl) appears in the lower left. Profiles have been averaged in time over roughly 200 days 
in each instance. A clear transition between solar-like and anti-solar rotation occurs between the middle and bottom rows (i.e. between 
the B and C series cases.) 


increasing along the horizontal axis. The upper row is 
comprised by simulations from Series A, and the lower 
row by Series C simulations. Convection models from Se¬ 
ries A and B each exhibit solar-like differential rotation, 
whereas series C tends to exhibit anti-solar rotation, with 
case C3 (lower right) possibly representing a transition 
point between the fast equator and fast pole regimes. 
Variation of the differential rotation with v is present, 
but the effect is greatly reduced with respect to that of 
k. This can be attributed to the v~ x k~ 2 scaling of the 


Rayleigh number noted above. 

As noted above, the key parameter in determining the 
transition between mean flow regimes is the Rossby num¬ 
ber Ro. Fig. [3] shows the relative differential rotation, 
AQ/Qq = (f } eq — Qpoie )/^Cb i n all simulations versus Ro, 
exhibiting a clear transition at R 0 ~ 0.17. This transi¬ 
tional value is lower than the R n y 1 value q uoted in 
previous papers (e.g. Gastine et al. 2013, 2014). We at¬ 
tribute this discrepancy to the ambiguity m defining Ro. 
Here we use the depth of the CZ as the length scale L 
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a 
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Ro 


Fig. 3.— Pole-to-equator differential rotation contrast AQ/Q as 
a function of Rossby number Ro for the A — C series of simulations. 
Transitional case B0.5 is indicated with a hollow diamond. The 
slightly more turbulent case B0, which is transitioning to an anti¬ 
solar state, is indicated by the red triangle. Ro was computed at 
mid-convection zone in all cases. Equatorial and polar rotation 
rates were calculated by averaging zonal velocity in both longitude 
and depth, and then over a latitude range of 0°N-20°N latitude 
and 70°N-90°N respectively. Anti-solar cases tend to arise when 
Ro is greater than about 0.163. 



and the rms value of the non-axisymmetric velocity in 


the mid-CZ v' 


as the velocity scale U. By contrast, 


(Gastine et al. 2014) use the volume-averaged value of 
v' rrns for U and a value of L based on the peak of the 
kinetic energy spectrum. This tends to shift U higher 
and L lower, resulting in a larger value of Ro. We refer 
the reader to ( Gastine et al.|2 014) for a discussion of how 
the transitional value of Ho depends on precisely how Ro 
is defined. 

The meridional circulation profiles corresponding to 
the simulations shown in Fig.[2]are shown in Fig. @ They 
exhibit a correlation with the differential rotation similar 
to that discussed in §4.1| Single-celled profiles of circula¬ 
tion appear in all of the C-series cases (i.e. those systems 
with fast poles), while those of the A-B series are multi¬ 
cellular in the equatorial regions. 

Convection in the solar and anti-solar regimes exhibits 
noticeably different morphologies as well. Fig. [5] il¬ 
lustrates the radial velocities in the upper and mid¬ 
convection zones of case Cl ( a,c ) and case A3 (b,d). Case 
Cl, which possess a anti-solar-like differential rotation, 
also exhibits velocity amplitudes that are roughly twice 
that of A3. Elongated convective cells in the equatorial 
regions tend to tilt in the negative ^-direction for case 
Cl. In addition to lower velocity amplitudes, Case A3 
possesses much more prominent columnar structuring of 
its convection, with a tilting evident in the positive in¬ 
direction. 

An additional effect is evident when comparing the 
flows of case A3 to those of case Cl, namely that the con¬ 
vection exhibits a much broader range of spatial scales. 
Case C3 is thus both more turbulent and more strongly 
driven. 


4.3. A Transitional Regime 


One might ask if the transition in mean-flow regimes 
is purely a result of rotational constraint, or if the level 
of turbulence might also play some role. In particular, 
by admitting much smaller spatial scales into the sim¬ 
ulation, are the correlations of v r and v$ necessary to 
establish a fast equator being disrupted? In other words, 
might the Reynolds number R e = ULjv play a role in 
addition to the Rayleigh number? In order to explore 
these questions we have extended the simulations of Fig. 
[2] to include the case B0.5. This case possesses a k identi¬ 
cal to the other B-series cases, but a v that is lower than 
any of of the C-series. Case B0.5 was initiated by lower¬ 
ing the diffusivity in the mature case Bl, which had itself 
been evolved for a total of two viscous diffusion times. 

An overview of the flows in Case B0.5 is shown in Fig. 
M The convective patterns are small-scale in nature, like 
C3, but do not exhibit the retrograde tilting character¬ 
istic of the patterns in this case. Case B0.5 possesses 
a decidedly solar-like differential rotation, but, surpris¬ 
ingly, seems to exhibit a single circulation cell within 
each hemisphere. Moreover, the flow amplitudes indi¬ 
cated near the surface are similar to those of case C3, 
as are the mid-convection zone velocities (not shown). 
Case B0.5 does possess a somewhat lower Rayleigh num¬ 
ber, and a different Prandtl number than case C3, how¬ 
ever, and so the differences between these two cases are 
likely to be subtle ones. It would seem that case B0.5 is 
just straddling the transition between the solar-like and 
a nti-solar regime s of differential rotation . 

Gastine et al.| ( 2014| ) and Kapyla et al. (2014) have 


recently shown that the behavior of the differential ro¬ 
tation near the solar/anti-solar transition can depend on 
the history of the simulation, representing a type of hys¬ 
teresis. For a given transitional i?o, the simulation may 
exhibit either solar or anti-solar DR. Though we have ob¬ 
served unusual mean flows near the transition with case 
B0.5, we have not searched explicitly for this bistability 
and hysteresis. It is possible that case B0.5 itself is slowly 
undergoing a transition from solar-like to anti-solar dif¬ 
ferential rotation. However, this case has been run for 
18,000 days, about a factor of four longer than the time 
scale for viscous diffusion across the convection zone, and 
we see no signs yet that it is undergoing such a transi¬ 
tion. The globally-averaged energy densities associated 
with the differential rotation (DRKE) and meridional cir¬ 
culation (MCKE) for case B0.5 are shown in Fig. frl The 
mean flows have varied by as much as 30% over this in¬ 
terval, but with no apparent long-term trends. We will 
discuss the significance of case B0.5 further in ^7| 


5. MAINTENANCE OF MEAN FLOWS 

In this section we argue that the distinct mean flow 
regimes described in 0 arise from a change in the nature 
of the convective angular momentum transport. Further¬ 
more, we argue that it is the convective angular momen¬ 
tum transport that largely determines the mean merid¬ 
ional flow profiles in simulations as well as stars. 


5.1. Dynamical Balances and Gyroscopic Pumping 

In an anelastic (low Mach number) system, the mass 
flux pv is divergenceless and the time evolution of the 
mean meridional flow p (v m ) is governed by 
vorticity equation, which may be written as 


t he zo nal 
(Miesch & 
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Contours^o^neridiona^irculatior^strearnline^withstreamftinction underlaVj^indicatedwithin 
each_^anei^TheJ^out and ti me averagin g are as i n Fig. [2] R ed ton es denot e cou nter-clock wise llow, an d b lue t o nes clockwise. As with 

cases where the rotation rate is varied, multi-cellular circulations arise whenever solar-like differential rotation protiles are present. 


HindmanpQll 


hereafter MH11) 


d_ 
d t 


= A 


on 2 

dz 


9 d(S) 
rC P d0 


(14) 


where = (Vxv m )»</> is the zonal component of the 
mean vorticity, A = r sin 0 and z = r cos 0 are the radial 
and axial coordinates in a (rotating) cylindrical coordi¬ 
nate system, and = Qq + A -1 (v^) is the total angu¬ 
lar velocity, including uniform and differential rotation 
components. We have adopted the convention that a 
subscript m on a vector denotes the radial r and latitu¬ 


dinal 0 components of a vector. We have focused here 
on the two dominant contributors to the meridional forc¬ 
ing, namely the inertia of the differential rotation and 
the baroclinicity of the mean stratification, which make 
up the first and second terms on the right-hand-side of 
(14) respectively. These are sufficient to illustrate the 
essential dynamics. We have neglected the meridional 
components of the Reynolds stress ((pv^v^); primes 
denote fluctuations about the longitudinal mean), the 
baroclinicity of thermal fluctuations, the viscous stress, 
and quadratic terms in v m . All of these are thought to 
be negligible in simulations and stars. We have also ne- 
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Fig. 5.— Radial velocity at one instant in time following equilibration for cases Cl (left column) and A3 (right column). The upper and 
lower rows correspond to radial velocity near the surface and in the mid CZ respectively. Convective cells in the solar-like case A3 exhibit 
a prograde tilt in the ^-direction. Those in anti-solar case Cl are tilted in the opposite sense. The limits of the color table normalization 
are indicated adjacent to each image, with units provided in m s -1 . 

cous stresses can ultimately determine the global struc¬ 
ture and amplitude of the flow. 

There is an alternative way to break the degeneracy of 
eq. (15) that provides a better explanation for how the 


glected the Lorentz force, which is less justified in stars 
but appropriate here since the simulations we consider 
are non- magnetic. Full expressions are given in MH11 
(see also Kit chat inov fe Rudiger | |1995 Rempel|2005 ). 

In a steady state eq. (|T1| gi ves the familiar thermal 
wind balance (TWB) equation ( Kitchatinov & Rudiger 

19951 Elliott et al.|2QQQ||Robinson & Chan 2 U(J 1 [|B " 


Toomre 20021 Rempel 


et al. 2009; Brun et at 


20051 Miesch et al. 2006[ Balbus 


2010 ) 




dtt 2 

dz 


9 d(S) 

A rC P 30 


(15) 


Here the two dominant meridional forcing terms bal¬ 
ance, and latitudinal thermal gradients can sustain non- 
cylindrical rotation profiles as in the Sun (891/dz ^ 0). 

It is important to note that in a steady state the merid¬ 
ional flow itself drops out of equation ( p~5] ) . Thus, it can¬ 
not be used to determine the MC profile even if other 
quantities such as 91 and S are known. This ceases to 
be the case if other terms are included in the balance, 
most notably the meridional components of the Reynolds 
stress. In mean field theory, this term is often represented 
as a turbulent diffusion. In short, inclusion of a turbu¬ 


lent diffusion term breaks the degeneracy of eq. (15) with 


respect to the meridional flow and can largely determine 
the MC profile that is achieved in many mean-field mod - 


els dKiiker et al.||2011| |Kitchatinov||2012| |Dikpati||2014|) . 
However, such solutions are sensitive to the nature of 
the imposed parameterizations, such as the depth de¬ 
pendence of the turbulent viscosity and departures from 
TWB in the upper and lower boundary layers, where vis- 


MC is established in our 3-D, time-dependent convec¬ 
tion simulations. The mechanism is known as gyroscopic 
pumping and can be illustrated by considering the zonal 
force balance which can be expressed as (MH11) 


p(v m ) -V£ = T * 


VF 


RS 


(16) 


where C = A 2 11 is the specific angular momentum, T is 
the net axial torque given by 

T = -V- [p\ (v>;> - -pv A 2 VO] = -V. [F *5 + F vd } , 

(17) 

and Frs is the convective transport of angular momen¬ 
tum, given by 


F RS = pX (w' m v'^) 


(18) 


In stars and in recent high-resolution convection simula¬ 
tions such as those presented here, the Reynolds stress 
component Frs generally dominates over the viscous 
component Fvd • The dynamical balances (15) and (|16|) 


should be regarded as temporal averages as well as lon¬ 
gitudinal averages, concerning the persistent mean flows 
that exist amid stochastic fluctuations. 

In convection simulations and in helioseismic inver¬ 
sions, V£ is oriented away from the rotation axis, parallel 
to the cylindrical radius VA (MH11). Thus, a conver¬ 
gence of the convective angular momentum flux (T > 0 ) 
will induce a meridional flow (v m ) directed away from 
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Fig. 6.— Sampling of flows for transitional case BO.5. (a) Radial 
velocity snapshot near the upper boundary, Upflows are indicated 
in yellow, downflows in blue. ( b ) Accompanying differential ro¬ 
tation (fi — Q©), time-averaged over 1800 days at the end of the 
simulation, with prograde rotation shown in red, and retrograde 
rotation in blue, (c) Meridional circulation streamlines for the 
same time period. Clockwise flow is indicated by blue underlay, 
and counter-clock wise flow by red underlay. This case demon¬ 
strates a prograde rotating equator in the absence of multi-cellular 
meridional circulations within each hemisphere. 


the rotation axis and a divergence ( J u < 0) will induce a 
flow toward the rotation axis. This is the phenomenon 


of gyroscopic pumping as discussed 

[ in a solar context by 

MH11 (see also Haynes et al. 1991 

Mclntyre| 1998[ 

2007 

Garaud & Arreguin 2009; Garaud < 

k Bodenheimer 

2010 

Miesch etal. 2012). 


nent of the rotational shear dVt/dz , but it can be sus¬ 
tained even if the steady-state rotation profile is strictly 
cylindrical; for an analytic illustration, see Appx. B of 
MH11. Indeed, this is also demonstrated in the simu¬ 
lations presented here; many exhibit Q profiles that are 
nearly cylindrical yet sustain persistent, well-established 
meridional flow profiles (e.g. Fig. [I]). Intuitively, one can 
attribute this to the efficiency ana robustness of Taylor- 
Proudman balance; in the barotropic limit (d (S) /d0 = 
0) any axial variation of the net torque {dJ-/dz ^ 0) 
will induce a transient axial shear (dfl/dz ^ 0) that 


will almost immediate!^] be wiped out by an induced 
meridional flow, via equation (14). In this way, a steady 


meridional flow can be maintained that continually re¬ 
plenishes angular momentum extracted or imparted by 
the net torque T by advecting angular momentum across 
local C isosurfaces as expressed by eq. (16). 

Eq. (16) predicts a direct link between the meridional 
flow prome and the angular momentum transport by the 
convective Reynolds stress. This link is strongest for the 
low-Rossby number limit in which the uniform rotation 
component of dominates C. In this limit we have 


* (A ’^ ) = 2 [ HX ' Z ' )dZ ' ’ 


(19) 


where zb = ( R 2 — A 2 ) 1 / 2 (MH11) and T is the cylindrical 
mass flux streamfunction, defined by 

<9T Id 

(p y x) = ^ . (20) 

Recall that A is the cylindrical radius (A = rsin#). Note 
that T is related to the mass flux streamfunction W as 
follows 


_ ld(W) 
r d0 


( 21 ) 


In summary, there are two competing mechanisms for 
breaking the degeneracy of (15) and establishing the 
steady-state MC; we will refer to them as the meridional 
Reynolds stress (RS) and gyroscopic pumping. Both rely 
on the RS tensor 7 Z^ = ~p{y r i v 1 ^. However, the first 
mechanism relies on the meridional components of 7Z 
(7 Z rr , TZor , 7^0 r , 1Z re , 7 Zqq, 7£</>#), whereas the second 
relies on the zonal components (7£ r 0, 7 Zq$, 7 Z^). In 
the next section we will demonstrate that it is this lat¬ 
ter mechanism, namely gyroscopic pumping, that deter¬ 
mines the MC profiles in our convection simulations. 

However, before proceeding, we first elaborate on these 
two mechanisms within the context of mean-field mod¬ 
els. In our convection simulations we capture all of the 
RS components explicitly so both of these mechanisms 
occur naturally. Both mechanisms may also be captured 
in mean-field models by parameterizing the RS based 
on convection simulations as presented here or on alter¬ 
native phenomenological arguments or turbulence mod¬ 
els. As mentioned above, the meridional RS mechanism 
may be captu red by solving the steady-state zonal vortic- 
ity equati on (Kiiker et al. 2011 Kitchatinov||2012 Dik- 
pati||2014|) . However, in order to capture the gyroscopic 
pumping mechanism it is essential to include the zonal 
momentum equation and to follo w the time-dependent 
approach to the steady state Rempel (20061). 

Though both mechanisms can and have been realized 
in mean-field models, the gyroscopic pumping mecha¬ 
nism is less sensitive to imposed mean-field parameteri- 
zations. There are two reasons for this. First, the merid¬ 
ional MC mechanism relies on second-order departures 
from the primary meridional force balance, which in the 
bulk of the convection zone is TWB, eq. (15). By con¬ 


trast, the zonal RS enters into the lowest-order zonal 
force balance, eq. (16). The second reason is that the 


3 on a timescale r ~ (2Q) 1 ~ P ro t/{ 47r), where P ro t is the 
rotation period. For the Sun this is about 2.2 days. 
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Time (10 Days) 


Fig. 7.— Evolution of mean-flow energy densities for transitional case BO.5. The energy in the differential rotation (DRKE) is depicted 
in black, and that of the meridional circulation (MCKE) in red. DRKE exhibits large (roughly 30%) variations with respect to the mean 
over this interval, but no long-term trends are evident. 


meridional MC mechanism is sensitive to the functional 
form of the RS parameterization, 7£(\k,0). Thus, im¬ 
plementing results from a convection simulation would 
require parameterizing the computed RS tensor. How¬ 
ever, in the gyroscopic pumping mechanism one must 
only know the axial torque T. This may in principle be 
decoupled from the mean flows without the need for any 
explicit parameterizations. The zonal components of the 
RS do likely depend on the amplitude of the differen¬ 
tial rotation, however, and a realistic mean-field model 
should take this into account. 


5.2. Convective Angular Momentum Transport 

In §5.1| we argued that the meridional flow profiles in 
our simulations are largely maintained by the convective 
angular momentum transport Frs through the mecha¬ 
nism of gyroscopic pumping. In this section we assess 
this argument and its implications by exploring the de¬ 
tailed nature of Frs and how it relates to the meridional 
flow profiles we find in the simulations. 

In our discussion, we find it instructive to introduce a 
Helmholtz decomposition of the convective angular mo¬ 
mentum transport Frs as follows: 

Frs = Vy + Vx . (22) 


Such a decomposition is valid for any arbitrary two- 
dimensional (axisymmetric) vector. Its utility is appar¬ 
ent when we see that only the first term contributes to 
the gyroscopic pumping equation (16). More generally, 
only the divergent component of the angular momentum 
flux can contribute to the time evolution of H, and thus 
the mean flow profiles that are ultimately achieved. One 
can derive the transport potential x from Frs by taking 
the divergence of (22) and solving a Poisson equation: 


V 2 x = V.F* 5 . (23) 

Figure [8] demonstrates the nature of the convective an¬ 
gular momentum transport in two representative cases, 
chosen to illustrate the low and high Rossby number 
regimes (upper and lower rows respectively). The trans¬ 
port potential x is illustrated in the left column (frames 
a, /). This is derived from the divergence of the convec¬ 
tive angular momentum transport shown in frames (5, g ) 
via eq. (23) and helps to elucidate its origin. 

In the slowly-rotating regime the angular momentum 
transport at mid latitudes is radially inward. This can 


be seen from the nearly horizontal orientation of the x 
contours in Fig. [8]f, increasing from top to bottom (ra¬ 
dially inward Vy). Since the radial component of Frs 
vanishes at the (impermeable) top and bottom bound¬ 
aries, this implies a divergence of the convective angular 
momentum transport in the upper CZ and a convergence 
in the lower CZ, as seen in the red and blue layers of Fig. 


0 ?- 


At low latitudes, Vy turns equatorward (Fig.l^f), pro 


ducing a convergence of Frs at the equator (Fig. |8p). 
Though this convergence of the convective angular mo¬ 
mentum flux at the equator tends to establish a solar- 
like differential rotation, it does not succeed in doing so. 
Rather, it is offset by the radially outward transport of 
angular momentum by the induced meridional flow (Fig. 
[8]/). This induced meridional flow, together with the 
radially inward angular momentum transport at higher 
latitudes is responsible for the anti-solar differential ro¬ 
tation profile th at is ultimately achieved. We discuss this 
issue further in §5.3| 

This scenario emphasizes the subtle nonlinear, non- 


local nature of the dynamical balances discussed in £5.1 


The convective angular momentum transport Frs is es¬ 
sential for establishing the differential rotation VH but 
it does not uniquely determine the resulting Q profile. 
Rather, Frs induces a meridional flow (mediated by the 
Coriolis force) that plays an essential role in establishing 
the differential rotation. 

The role of the convective angular momentum trans¬ 
port in establishing the meridional circulation is illus¬ 
trated in Fig. [8^. As expressed by eq. (fl6|, a convergence 
of F rs (blue) induces a flow away from the rotation axis 
(white arrows) and a divergence of Frs (red) induces a 
flow toward the rotation axis (black arrows). Mass con¬ 
servation then requires these flows to form closed circula¬ 
tion cells, as reflected by the actual meridional flow pat¬ 
tern shown in Fig. |8^. Due to the nature of F#c( which 
can ultimately be traced to the x profile in Fig.j^f), the 
induced meridional flow has a single-celled nature, dom¬ 
inated by one large cell in each hemisphere that extends 
from the equator to a latitude of at least 70° and spans 
the entire CZ. These cells exhibit poleward flow in the 
upper CZ and equatorward flow in the lower CZ. The ad- 
vection of angular momentum by this induced flow nearly 
balances the convective angular momentum transport as 
expressed by eq. (16). This can be seen by comparing 
Figs. Hb and [8]/. small imbalances are due to viscous 
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Fig. 8. — Maintenance of meridional flow. Illustrative examples are shown for the rapidly-rotating regime (upper row, case Cl) and the 
slowly-rotating regime (lower row, case A3), (a, /) Transport potential x> increasing from -10 26 g cm s -2 (blue) to 10 26 g cm s -2 (red). 
(6, g ) Divergence of the convective angular momentum transport, ranging from -10 7 g cm -1 s -2 (blue) to 10 7 g cm -1 s -2 (red), (c, h) 
same as frames ( b , g) but with arrows indicating the direction of the induced meridional flows. Slanted lines in the northern hemisphere 
delineate high and low-latitude regimes, (d, i) Meridional flow shown as streamlines of the mass flux T. Red denotes clockwise flow ('F > 0) 
and blue denotes counter-clockwise flow ('F < 0). Values of \F range from ±4 x 10 11 g cm -1 s -1 in (d) and ±9 X 10 10 g cm -1 s -1 in (i). 
(e, j) Angular mom ent um transport by the meridional flow which approximately balances the Reynolds stress divergence in frames ( b , g), 
as expressed in eq. ( |16| ). Color table as in frames (6, g). 


torques [eq. ©] and unsteady fluctuations. 

Similar arguments hold for the rapidly-rotating regime 
in the upper row of Fig. (|8|. Here the x contours at high 
latitudes are again nearly norizontal, signifying a radially 
inward angular momentum transport. However, outside 
the tangent cylinder, the angular momentum transport is 
cylindrically outward, away from the rotation axis. This 
is reflected by the cylindrical nature of the x contours 
at low latitudes. There is also an equatorward contribu¬ 
tion to Vx, particularly at mid-latitudes near the outer 
boundary. This x profile produces a Reynolds stress di¬ 
vergence pattern that reverses sense at high and low lat¬ 
itudes, as seen in the red and blue regions of Figj8^. At 
high latitudes the pattern is similar to the high-Ro case, 
with a divergence and convergence in the upper and lower 
CZ respectively. At low latitudes this reverses, exhibit¬ 
ing a convergence in the upper CZ and a divergence in 
the lower. 

The implications of this x profile for the meridional 
circulation are profound. Inside the tangent cylinder, a 
single cell is established much like the high-Ro case (Figs. 


Hb i). However, outside the tangent cylinder, a multi¬ 
cell profile is induced with 2-3 distinct cells spanning 
the CZ (Fig. [§i) . This can be attributed largely to the 
cylindrically-outward angular momentum transport near 
the equator which induces upward and downward flow 
in the upper and lower CZ respectively (Figs. [8]c). As 
in the high-Ro regime, the advective angular momentum 
transport by this induced meridional flow largely bal¬ 
ances the convective angular momentum transport (Figs. 
[8b, e), with small departures primarily due to the contri- 
oution of viscous diffusion to the zonal force balance. 

The dramatic difference between the convective angu¬ 
lar momentum transport at high and low latitudes can 
be appreciated by considering the structure of the con¬ 
vective flow and how it is influenced by rotation. At high 
latitudes, convection is dominated by a quasi-isotropic, 
interconnected network of downflow lanes near the sur¬ 
face that break up into more isolated lanes and plumes 
deeper down (Fig.j5j). As these downflows travel down¬ 
ward, the Coriolis force deflects them in a prograde direc¬ 
tion, inducing a negative correlation between v' r and v f . 
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that produces an inward angular momentum transport 
(F r S ). 

At low latitudes, outside the tangent cylinder, the pre¬ 
ferred convection modes for low Ro are banana cells , 
columnar convective rolls that are alig ned with the rota¬ 
tion axis near the equator (e.g., Miesch|2 QQ5~). In simula¬ 
tions with large density stratification and moderate val¬ 
ues of Ro, these banana cells are manifested as a promi¬ 
nent north-south alignment of downflow lanes near the 
equator that often trail off eastward at higher latitudes 
where they are distorted by the differential rotation (Fig. 
|5^; see also Miesch et al. 2008). 

In their simplest manifestation at low Ro and moder¬ 
ate density stratification, banana cells are approximately 
aligned with the rotation axis with little variation in the 
axial (z) direction. The combined action of the spherical 
geometry, the density stratification, and positive nonlin¬ 
ear feedback from the differential rotation they establish, 
tends to produce a systematic tilt such that cylindri- 
cally outward flows (v\ > 0, where A is the cylindrical 
radius) are deflected west ward and cylindrically inward 
flows (v\ < 0) eastward (|Busse|2002||Aurnou et al.|2007|). 


This establishes a negative ( v' x v correlation that is re¬ 
sponsible for the cylindrical orientation of the x contours 
at low latitudes in Fig.[8]a (i.e. the cylindrically outward 
angular momentum transport). The effects of density 
stratification, the outer spherical boundary, and moder¬ 
ate (but still low) Ro tend to cause banana cells to bend 
away from the z axis toward a more horizontal orienta¬ 
tion, with axes more parallel the 0 dimension. This tends 

to establish positive (vgV^ correlations (in the northern 

hemisphere, negative in the so uthern) that transport an¬ 
gular momentum e quator ward (Gilman 1983 Glatzmaier 
1984l|Miesch|[2QQ5|). 


Thus, banana cells can account for both the cylindri¬ 
cally outward angular momentum transport at low lat¬ 
itudes in the low -Ro regime (Fig. If) and the equator- 
ward angular momentum transport at low latitudes in 
the high -Ro regime (Fig. [8]f). Though the region out¬ 
side the tangent cylinder contains the banana cells, a 
close look at Fig. [8|a indicates that this region does not 
strictly delineate the transition between the high and 
low-latitude regimes for the angular momentum trans¬ 
port. Rather, the region of inward transport at high 
latitudes shifts equatorward with increasing Ro, as in¬ 
dicated schematically with the solid black lines at mid 
latitudes in Figs. [8]c and[8^. 

This equatorward shift of the transition between high 
and low-latitude behavior can be understood by consid¬ 
ering a downflow plume spawned from the upper bound¬ 
ary layer. Though such a plume is subject to a range of 
forces including pressure gradients, buoyancy, nonlinear 
advection and viscous diffusion, it is instructive to con¬ 
sider a ballistic trajectory subject only to the Coriolis 
force. Then a downflow plume with initial velocity — v p r 
will be deflected in a prograde direction with an initial 
radius of curvature given by 


(24) 


dv^/dt 


2f2o sin 9 


measure of the local, latitude-dependent Rossby number 
experienced by a downflow plume: 


p = l£ _ V _P 

8 D 2QqD sin 9 


(25) 


Intuitively, eq. (25) can be interpreted based on whether 
or not a downflow plume on a ballistic trajectory will tra¬ 
verse the CZ before being deflected by the Coriolis force. 
If Rq » 1 then the plume could in principle (in the ab¬ 
sence of other forces) reach the base of the CZ with only 
a small prograde (positive </>) deflection. This Coriolis- 
induced deflection reflects the tendency for the plume 
to conserve its angular momentum conservation and it 

is responsible for the negative correlation that 

transports angular momentum inward at high latitudes 
as seen in Figs. [8]a and /. 

For Rq « 1 the plume will not make it to the base of 
the CZ (again, considering the simplified case of a bal¬ 
listic trajectory). Rather, the convective flow will have 
to restructure itself in order to provide the requisite heat 
transport. This is the realm of banana cells. This sug¬ 
gests that the transition between the two regimes should 
occurs where Rq ~ 1. However, we found above that 
the actual transition occurs at a lower threshold Rossby 
number of R t ~ 0.16 when Ro is defined based on the 
rms velocity and the depth of the CZ (cf. Fig. [3|. In 
the present context this implies that the transition to 
the rapid-rotation regime occurs when a ballistic plume 
is deflected well before it reaches the base of the CZ; 
r c < RfD. This suggests a transition colatitude of 


0o 


2 q 0 d 


Ro 

~Rt 


(26) 


where Ro is the global Rossby number defined in ^l]and 
where we have assumed that v v ^ v rms . For the solar- 
like case A3 (Ro = 0.11), this corresponds to a transition 
latitude of 50°. This estimate agrees well with the actual 
transition latitude indicated in Fig. 

Equation (26) should be regarded only as a loose rule 


of thumb for accounting for why the transition between 
polar and equatorial regimes for the convective angular 
momentum transport shifts equatorward with increas¬ 
ing Ro. It formally breaks down for Ro > Rt where 
Re > Rt at all latitudes. Yet, even in this regime, ba¬ 
nana cells contribute equatorward angular momentum 
transport near the equator as seen in Fig. [8]f. 

In summary, we find that the angular momentum 
transport by the convective Reynolds stress, plays 

a central role in establishing the meridional flow profiles 
in our simulations. In particular, we find a transition 
from single-celled to multi-celled meridional circulation 
profiles in the high and low Ro regimes that is directly 
linked to a change in the nature of ¥rs- 

This conclusion is su pporte d by the recent high- 
resolution simulations of Hotta et al. (2014). There they 


find that inward angular momentum transport by small- 
scale convection in the surface layers induces a poleward 
meridional flow as envisioned by MH11. As in our sim¬ 
ulations, this is a clear demonstration of M C pr ofiles es¬ 
tablished by gyroscopic pumping (see also [7.2). 


The value of r c relative to R, the depth of the CZ, is a 


5.3. The Solar-Anti-solar Transition 
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In §5.2| we demonstrated that the convective angular 
momentum transport F^s, plays a central role in regu¬ 
lating the MC profile in both the high and the low Ro 
regimes. In this section we demonstrate that Frs also 
regulates the transition between the two regimes and that 
the induced MC plays a key role in this transition. 

We focus primarily on the anti-solar (high Ro) regime 
where the convective angular momentum transport is ra¬ 
dially inward at most latitudes (Fig.pV). As discussed in 
§5. 2 1 and as shown in Fig. [8]i, this induces a single-celled 
MC that is counter-clockwise in the northern hemisphere 
(NH) and clockwise in the southern hemisphere (SH). 
This MC transports angular momentum that offsets the 
convective angular momentum transport (Fig. [8]/). 

However, the MC induced by gyroscopic pumping also 
transports entropy. Since the stratification of the CZ is 
superadiabatic ( dS/dr < 0), the upflows and downflows 
at low and high latitudes tend to establish an equator- 
ward entropy gradient ( dS/dO positive in the NH and 
negative in the SH). This gives rise to a baroclinic forcing 
that further enhances the MC. Baroclinicity thus pro¬ 
vides a positive feedback that amplifies the MC estab¬ 
lished by gyroscopic pumping. Essentially, the mechani¬ 
cal forcing due to Frs triggers an axisymmetric convec¬ 
tion instability that feeds back on the mean flows by re¬ 
distributing both angular momentum and entropy. This 
axisymmetric mode is superposed on the preferred non- 
axisymme tric modes for convectio n in rotating spheri¬ 
cal shells (Miesch & Toomre 2009 e.g.). Though rota¬ 
tion plays an essential role in exciting this axisymmetric 
convective mode through the Coriolis-induced Reynolds 
stress, it may have some bearing on similar large-scale 
circulations observed in laboratory experiments and sim¬ 
ulations of n on-rotatin g turbulent Rayleigh-Benard con¬ 
vection (A hlers et al.||20 09). 

This is demonstrated in Fig. [9| which illustrates how 
the anti-solar differential rotation in Case Cl is estab¬ 
lished. Shown in blue is the time evolution of L, which is 
the angular momentum of the north polar region relative 
to the rotating reference frame: 


Lp — 


rr 2 r 6 p r 2n 

f f f P r3 

Jn Jo Jo 


pr 6 sin Ov^drdOdcj) 


(27) 


where 0q = tt/6, corresponding to a latitude of 60°. In 
Fig- HO L is normalized by L 0 , the total amount trans¬ 
port poleward by meridional circulations over this inter¬ 
val. The stress-free boundaries exert no torques so the 
total angular momentum [obtained by replacing 0q with 
7r in eq. (27)] is conserved. Thus, the time evolution re¬ 
sults frornthe integrated angular momentum flux across 
a latitude of 60°. This flux includes contributions from 
the convective Reynolds stress, the meridional circula¬ 
tion, and the viscous diffusion, represented in Fig. ibya 
dashed line, a solid black line, and a red line respectively. 

The first thing to note in Fig. [9] is that the convec¬ 
tive angular momentum transport is not strictly radial. 
Rather, it includes a weak equatorward component that 
acts to spin up the equator even at mid latitudes. Over 
the first 4000 days of evolution, this is approximately bal¬ 
anced by the meridional circulation, which transports an¬ 
gular momentum poleward. However, eventually the MC 
overwhelms the Reynolds stress, spinning up the poles 
until viscous diffusion sets in to help limit the growth of 



Fig. 9.— Polar-spinup in anti-solar case Cl. Total angular mo¬ 
mentum transported poleward across 60° N latitude, as a function 
of time, is shown in blue. The individual contributions to this 
transport by the meridional circulations (black line), convective 
Reynolds Stresses (dashed line), and viscous stresses (red line) are 
indicated as well. Polar spin-up in case Cl arises from an im¬ 
balance in transport by Reynolds stresses, which work to spin up 
the equator, and transport by meridional circulations, which work 
to spin up the poles. Equilibration occurs around day 6000 as vis¬ 
cous stresses associated with strong rotational shear become strong 
enough to help counter the effect of meridional circulation. 


the cyclonic polar vortex. 

Thus, it is the meridional circulation that ultimately 
establishes the anti-solar differential rotation, not the 
Reynolds stress directly, although the latter induces the 
former. The transition from solar to anti-solar differen¬ 
tial rotation is triggered by the inward angular momen¬ 
tum transport but mediated by the strong, single-celled 
MC. This stems from the tendency for the MC to homog¬ 
enize angular momentum, spinning up the poles relative 
to the equator (Paper 1) and is consistent with the pole- 
ward angular momentum transport by the MC re porte d 


in many global convection simulations (IBrnn & Toomre 
2002 Miesch et al.||2008 Kapyla et al.||201lT 
The role of the 1VIC in establishing anti-sol< 


establishing anti-solar DR pro¬ 
files h as previou sly bee n emph asized in mean-field mod¬ 
els by Kitchatinov & Rudiger (2004|) (see also Rudiger & 
Kitchatinov 2007). However, in these models the 1VIC is 
not established by convective transport. Rather, other 
processes such as the suppression of high-latitude con¬ 
vection by polar starspots or tidal forcing from a bi¬ 
nary companion induce baroclinic torques that in turn 
induce global circulations. By contrast, in our convec¬ 
tion models, the MC is established by a radially inward 
convective angular momentum transport and reinforced 
self-consistently by baroclinicity. 


6. CONVECTIVE HEAT TRANSPORT IN SOLAR AND 
ANTI-SOLAR CASES 


In §5.If§5.2| we argued that it is the convective an¬ 
gular momentum transport that ultimately determines 
the gross differential rotation contrast AD as well as the 
MC profile. However, we have also seen that baroclin¬ 
icity can play an important role as well, both in estab¬ 
lishing anti-solar DR profiles by spinning up the poles 
(j|5.3|) and in shaping the D contours through thermal 
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Fig. 10. — Energy flux balance in cases B1 (a; our most turbulent 
solar-like case) and Cl (6; the most strongly driven anti-solar case). 
Fluxes have been integrated over the shell to yield a luminosity. 
Enthalpy flux is shown in blue, radiative flux in red, conductive 
flux in solid black, and inward kinetic energy flux in dotted black. 
A dashed line has been plotted at unity for reference. Convection 
tends to transport additional heat in the anti-solar case relative 
to the solar case, compensated for by an increased inward kinetic 
energy flux. 


wind balance, eq. (15). Though most convection and 


mean-field models agree that TWB is resp onsible for the 
conical Cl profiles in the solar CZ (j ]5.l| ), the origin of 
this baroclinicity is not well understood. Early models 
attributed it to the influence of rotation on convective 


heat transport ([K itchatin ov fc Riidig er|| 19951 |Robinson 
& Chan||2001| |Brun &; Toomre||2002p but thermal cou¬ 
pling to the t achocli ne may also p lay an important role 
(Rempel|2005 Miesch et al.j2006 ). Here we consider the 
former by investigating the inhomogeneous nature of the 
convective heat flux in the solar and anti-solar regimes. 

A sense of the energy flux balance in the rapidly- 
rotating and slowly-rotating regime may be found by 
examining Figure [K)| Plotted in Figure [To| are the con¬ 
vective enthalpy flux F e , the radiative heat flux F r , the 
kinetic energy flux iA, and the conductive heat flux F c , 
defined as 

F e = pc p v r T\ (28) 


„ _ - dT 

r r — K r pC p ^ , 

1 


F k = ^pv r v 2 , 


(29) 

(30) 


Fig. 11.— Variation of entropy in the (a) low Rossby number 
(case A3) and ( b ) high Rossby number (case Cl) regimes. The 
spherically symmetric component of the entropy perturbations is 
plotted in black for each case. Cuts at the equator and a latitude 
of 85° have been plotted in blue and red respectively. Variations 
about the spherically symmetric mean have been amplified by a 
factor of 3 for visibility. Convection is more efficient at high lat¬ 
itudes in the low Rossby number regime, signified by the steeper 
entropy gradient near the surface and the nearly adiabatic inte¬ 
rior in those regions. For the low Rossby number regime of (6), 
convection is more efficient at low latitudes. 


and 




(31) 


Each of these fluxes has been integrated over horizon¬ 
tal surfaces and normalized by the solar luminosity. The 
outward convective transport of heat and inward trans¬ 
port of kinetic energy dominate the flux balance at mid¬ 
convection zone in either system. The rapidly rotating 
system B1 (Figure [l0|a) has somewhat less efficient con¬ 
vection than the slowly-rotating system (Figure [lOp), 
and thus possesses a stronger contribution from conduc¬ 
tion at mid-convection zone. This behavior is typical 
of other rapidly rotating systems; namely a strong ro¬ 
tational constraint tends to inhibit the convective effi¬ 
ciency. The latitudinal averaging that has gone into Fig- 


10 belies significant differences in the convective state 


achieved by these systems, however. 

The convective transport of heat achieved in either the 
solar or anti-solar regimes is highly anisotropic latitudi- 
nally, with the nature of that anisotropy depending on 
which rotational state the convection is in. Due to the 
impenetrable outer boundary present in these systems, 
transport of heat across the upper boundary is accom- 
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spherically symmetric mean subtracted at each radius, for solar-like 
case A3 and anti-solar case Cl. The solar-like cases tend to develop 
polar regions that warm with respect to the equator. By contrast, 
the anti-solar case exhibits a warm equator. These longitudinal 
averages have also been time-averaged over a period of roughly 
200 days. 


plished by conduction. It is in this region where the 
super-adiabaticity of the system is established, resulting 
from the development of a conductive boundary layer. 
As we have chosen to fix the value of entropy at the up¬ 
per boundary, the entropy gradient there is free to vary 
in time, and also in latitude. Fig. [Tl] illustrates the typ¬ 
ical thermal background realized in the rapidly rotating 
(Fig. 11a; case A3) and slowly-rotating (Fig. [TT) fr; case 
Cl) regimes. The spherically symmetric entropy profile 
is plotted in black alongside cuts of the axisymmetric 
entropy profiles at the equator (red) and 85° (blue). At 
high latitudes, super-adiabaticity in the entropy gradi¬ 
ent is confined to a thin boundary layer for the solar-like 
case, implying efficient convection. The situation is the 
opposite at low latitudes, where much of the convection 
zone is super-adiabatic, albeit with a lower overall adia- 
bat. This situation, with efficient polar convection and 
inefficient equatorial convection is typical of the rapidly- 
rotating regime and arises from the prominent role played 
by the Coriolis force, which tends to inhibit radial down- 
drafts more efficiently at the equator than at the poles. 

The slowly-rotating regime exhibits an even more pro¬ 
nounced dichotomy between pole and equator, but in¬ 
verted. Here convection is most efficient at the equa¬ 
tor, with a super-adiabatic stratification at high latitudes 
(Fig. [TT)fr). Here the convection is suppressed by the 
strong prograde polar vortex that is characteristic of the 
anti-solar regime. Down-welling plumes are quickly de¬ 
flected in the prograde direction by the strong rotational 
shear and the cylindrically outward Coriolis force. The 
polar vortex also acts as a transport barrier for merid¬ 
ional momentum, suppressing the meridional flow as can 
be seen by examining the C cases in Fig. [4j 

The thermal anisotropies realized in the two different 


regimes have implications for the establishment of merid¬ 
ional circulation as well. Profiles of the axisymmetric en¬ 
tropy perturbations are plotted in Fig. [12] where, in order 
to emphasize the latitudinal variation oFthe entropy per¬ 
turbations, we have subtracted out the spherically sym¬ 
metric mean. In case A3 (Fig. [l2]a), the polar regions 
are warm with respect to the equator. In isolation, this 
effect alone would drive a circulation that is equatorward 
at the surface. This is not the nature of the circula¬ 
tion realized in the rapidly rotating regime. In these 
simulations, the warm poles thus serves to suppress the 
meridional circulations realized. The meridional circula¬ 
tions in the rapidly-rotating regime are thus not driven 
thermally. The same can be said of the Sun (Paper 1). 

The situation at the surface of Cl is reversed, how¬ 
ever, and the effect is to augment any poleward circu¬ 
lations already established in the near-surface layers by 
gyroscopic pumping. This reflects the excita tion of an 
axisymmetric convective mode as discussed in §5.3| Here 
baroclinicity is unable to achieve thermal wind balance 
and viscous forces are needed to establish equilibrium. In 
actual stars where viscous forces are negligible, this role 
may be played by magnetic fields, which would inhibit 
the development of a polar vortex. 


7. APPLICATION TO THE SUN AND STARS 
7.1. What Regime is the Sun in? 

What do our results imply about the meridional cir¬ 
culation in the Sun? Is it single-celled or multi-celled? 
Since its rotation profile is solar-like, by definition (All > 
0), we might readily place it in the rapidly-rotating 
regime and we might therefore expect that it should have 
a multi-celled profile. However, the situation is more 
subtle than this simple categorization suggests. The dif¬ 
ficulty lies with the conical nature of the Sun’s D profile. 

There is a clear difference between the “solar-like” dif¬ 
ferential rotation profiles in the low-Rossby number sim¬ 
ulations presented in ®] and the actual internal rotation 
profile of the Sun inferred from helioseismology. Like the 
Sun, the simulations in the rapidly rotating regime ex¬ 
hibit a fast equator and slow pole (All > 0). However, 
unlike the Sun, the simulations exhibit rotation profiles 
that are approximately cylindrical, such that Vfi is di¬ 
rected away from the rotation axis. By contrast, helio- 
seismic rotational inversions suggest that rotation pro¬ 
files are more conical in nature, such that Vfl is pri¬ 
marily equatorward, with nearly radial D is osurfaces at 
mid-latitudes as shown in Fig. |l3|a (see also Thompson 
et al. 2003 Howe] 2009). The presence of these conical 
profiles may signify that the Sun is near the transition 
between the two mean flow regimes. 

To illustrate why this is the case, we must first un¬ 
derstand why the D profile in the Sun is coni cal rather 
than cylindrical. As discussed in section 5.1, a diverse 
variety of modeling efforts including convection simula¬ 
tions, mean-field models, and theoretical arguments have 
converged on a common explanation. They attribute 
the conical nature of the solar rotation profile to ther¬ 
mal gradients, which contribute a baroclinic forcing to 
the zonal vorticity equation that can effectively offset 
meridional flows induced by the inertia of the differential 
rotation that would otherwise establish cylindrical pro¬ 
files in accordance with the Taylor-Proudman theorem. 
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This robust balance between baroclinicity and inertia is 
expressed by the thermal wind equation (15), which im¬ 
plies 


9H A „ 
0—- oc AS 
oz 


(32) 


where AS is the entropy difference between pole and 
equator. How the differential rotation profile changes 
as the global rotation rate 0 is increased ultimately de¬ 
pends on the scaling of AS with rotatio n rate. 

The mean-field models of Hotta & Yokoyama (2011), 
which include thermal coupling to the - oachocline, suggest 
that AS does not scale steeply enough with O to main¬ 
tain a conical differential rotation for fast rotators. These 
models exhibit an increasingly cylindrical differential ro¬ 
tation profile as D is increas ed. Similarly, simulations of 
G-type sta rs by|Brown et ah1(|2 008|) and recent mean-field 
models by Kiiker et al. [funjr also show nearly cylindri¬ 
cal alignmen t at rapicfrotation rates . However, in mean- 
field models ([ Kiiker & Rudiger 2005) and 3-D simulations 
(Augustson et al.|2012|) of convection in F-type stars with 
thinner convection zones, thermal gradients are able to 
scale with rotation gradients, maintaining a significant 
conical orientation over nearly an order of magnitude in 

a 


If D profiles do become more cylindrical at faster ro¬ 
tation rates, then the conical D profile of the Sun sug¬ 
gests that it is in a low-Rossby number (rapidly-rotating) 
regime but perhaps only marginally. What does this im¬ 
ply about the meridional flow profile? 

In section [5] we argued that the multi-cell merid¬ 
ional flow profile at low latitudes in rapid rotators can 
be attributed to cylindrically outward angular momen¬ 
tum transport by columnar convective modes (“banana 
cells”). This tends to produce not only a multi-cell cir¬ 
culation profile but also a cylindrical D profile with a 
prominent dTl/dr at the equator (see Fig. Tj) . 

The solar rotation profile exhibits a weax 90 /dr near 
the equator with the same sense (positive) as in the 
rapidly-rotating solar models. This may indicate that 
some cylindrically outward angular momentum transport 
is indeed occurring. However, the angular velocity con¬ 
trast AO across the CZ at the equator (10-20 nHz) is 
small compared to the contrast between the equator and 
mid-latitudes (50-60 nHz), suggesting that the convective 
angular momentum transport may be more equatorward 
rather than cylindrically outward. 

We do indeed find equatorward angular momentum 
transport by convection at low latitudes in the slowly 
rotating regime, as demonstrated by the nearly radial x 
contours adjacent to the equator in Fig. [§]f. This conver¬ 
gence of angular momentum flux at the equator tends to 
accelerate the rotation nearly uniformly throughout the 
CZ. Furthermore, through the mechanism of gyroscopic 
pumping, it establishes a single-cell meridional circula¬ 
tion profile, with radially outward flow at the equator 
that veers poleward at the surface. 

In the slowly rotating regime, the convergence of the 
convective angular momentum flux at the equator is not 
sufficient to establish a solar-like differential rotation pro¬ 
file. Rather, the radially inward angular momentum 
transport at higher latitudes and the meridional flow 
it establishes efficiently spin up the poles, producing an 
anti-solar profile. However, there might be situations in 


which equatorward angular momentum transport (as op¬ 
posed to cylindrically outward) can sustain a solar-like 
O profile. 

One such exa mple is our transitional case B0.5 dis¬ 
cussed in section |4.3| As shown in Figure [6j this case 
exhibits a solar-like O profile (AD > 0) with a single¬ 
cell meridional flow. The convective angular momentum 
transport in this case is shown in Fig. [l4j As in the 
anti-solar regime (Fig. |8]f), it is radially inward at high 
latitudes and equatorward at low latitudes. Yet here the 
O profile is solar-like. 

Another way to approach t his issue is to turn the prob¬ 
lem around, as described by Miesch (2005). Using the 
solar rotation profile inferred from nelioseismology, we 
can ask what Reynolds stress would be needed to main¬ 
tain it for a given meridional flow profile. Fig. [T3| shows 
the result for a single-celled profile. In other words, Fig. 


13 d shows the Reynolds stress that would be needed to 


establish a single-celled MC profile in the Sun, expressed 
in terms of the transport potential y. This looks notably 
similar to the x contours in the transitional case B0.5, 
shown in Fig.[l4|a. The latter are more radial at low lati¬ 
tudes but this is sensitive to the precise form of the MC; 
what matters here is the qualitative similarity. Also, the 
O profile in case B0.5 is not completely analogous to the 
Sun in that it is cylindrical rather than conical. How¬ 
ever, this discrepancy may have to do with the absence 


of a tachocline, as discussed in £7.3 


Our results suggest that the Sun may be near the tran¬ 
sition between the fast and slow rotating regimes. This 
is consistent with previous estimates of the Rossby num¬ 
ber of the Sun, which he just below unity (Paper 1). We 
predict that G-type stars rotating faster than the Sun 
should have cylindrical O profiles with multi-cell merid¬ 
ional circulation profiles and those rotating much slower 
than the Sun should have anti-solar O profiles and single¬ 
cell meridional flow profiles. However, the Sun is close 
enough to the transition that we cannot definitively say 
from these simulations alone whether it has a multi-cell 
or single-cell profile. 

It is interesting to note that if the mean flows of the 
Sun are in a transitional regime and thus atypical relative 
to other stars with faster or slower rotation rates, then 
this should be reflected in its observed patterns of mag¬ 
netic activity. Both the differential and the meridional 
circulation play an essential role in mos t rece n t dyn amo 
models of the solar cycle (e.g. 
these are atypical, then the so 


Charbonneau 2010). If 
ar cycle should be too. 
Although data on stellar cycles is limited and sometimes 
ambiguous, there is some evidence that the period of the 
solar cycle relative to its rotation rate is indeed atypical, 
lying betwee n the active and inac tive dynamo branches 
described by Bohm-Vitense (2007). 

Further insight must come from continuing solar and 
stellar observations, helioseismic inversions, and from the 
consideration of factors such as boundary conditions and 
magnetism, which can influence the solar/anti-solar tran¬ 
sition and the structure of mean flows. We address these 
issues in the following subsections. 


7.2. The Near-Surface Shear Layer 

The simulations reported on here lack a notable com¬ 
ponent of the solar convection zone: a near-surface shear 
layer (NSSL). This region of shear, located within the 
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Fig. 13.— Inferred angular momentum transport for the Sun, assuming a single-cell meridional circulation profile (after IMiesch|| 2005l>. 
(a) Angular velocity Q based on RLS inversions of GONG data for four non-overlapping intervals in 1996, provided by R. 11 owe (|Howe 
|et al. 2000 Sch ou et al.|2 002). (6) Hypothetical single-cel led me ridional circulation profile, shown in terms of the streamlines of the mass 
flux, used in the flux-transport dynamo models of|Dikpati (2011). Blue denotes counter-cl ockw ise circulation; the amplitude is arbitrary for 
our purposes, (c) Mean axial torque required to sustain these mean flows, from equation ( |16| ): T = ~p (v m ) •'VC. Red denotes a divergence 
of the angular momentum flux (J 7 < 0) and blue denotes convergence, (d) Corresponding x contours, as in Fig. [8] 



Fig. 14.— Convective angular momentum transport in the tran¬ 
sitional case, B0.5. Shown are (a) the transport potential x and ( b ) 
the negative divergence of (cf. Fig. (8] frames a, b , /, and g). 
Color tables are as in Fig.[8]with saturation levels of (a) 1.4xl0 26 
g cm s -2 and ( b ) 2.1 xl(F g cm -1 s -2 . 


outer 5% of the solar interior, is characterized by a de¬ 
crease in an gular velocity with radius (dfl/dr < 0) of 
about 2-4% (Thompson et al. 2003 Howe] 2009). Our 
simulations lack an IN SSL comparable to the bun likely 
because of the simplified upper boundary condition and 
the lack of small-scale surface convection (recall from ^3] 
that our simulations only extend to 0.965R). However, 
the insight gained from these simulations as discussed 
in HI can help guide the interpretation of helioseismic 
inversions. 

As discussed by MH11, the existence of the solar NSSL 
is often attributed to the inward transport of angular 
momentum by small-scale convection in the solar surface 
layers where the effective Rossby number is high. The 
phenomenon of gyroscopic pumping links the differential 
rotation and meridional flow in the NSSL and can, at 
least in part, account for the poleward flow near the so¬ 


lar surface. In particular, the inward angular momentum 
transport by solar surface convection gives a retrograde 
torque (flux diver genc e) that sustains the poleward flow 
through equation ( |16| ) . Much of this inward angular mo¬ 
mentum transport likely arises from convective plumes 
that are driven by photospheric cooling and that are de¬ 
flected in a retrograde direction by the Cor iolis force as 
they propagate downward, as described in §5.2| If these 
plumes penetrate to a radius of r p , then one would ex¬ 
pect a requisite convergence of the angular momentum 
flux which would then induce an equatorward flow near 
r p as expressed by the gyroscopic pumping equation (16). 

Thus, the recent detection of a shallow equatorward 
return flow at r ^0.91-0.9377 may reflect the penetration 
depth of convective plumes that originate in the pho¬ 
tosphere. This is consistent with the interpretation of 
Hathaway (2012b) who attribute the shallow return flow 
to the presence of the NSSL and the Reynolds stress asso¬ 
ciated with supergranular convection. Furthermore, this 
may account for w hy th e shallow return flo ws reported by 
Hathaway (2012b) and Zhao et al. (2013) reverse direc- 
tion along surfaces that are nearly horizontal while the 
reversal surfaces in our multi-celled low-Rossby number 
simulations are more cylindrical (see Fig. [l]f). It will be 
interesting to see if future inversions show any signs of 
cylindrical alignment in deeper layers where the Rossby 
number is low. 

More generally, helioseismic inversions that focus on 
the NSSL region may reflect local dynamics that have lit¬ 
tle bearing on the deeper dynamics responsible for main¬ 
taining the solar-like differential rotation (Afi > 0) and, 
by extension, the deep meridional flow. Since the process 
of gyroscopic pumping is intrinsically nonlocal (MH11), 
the presence of the NSSL may influence the meridional 
circulation profile throughout the convection zone. How¬ 
ever, the lower CZ is certainly not passive; this much is 
clear from the nature of the H profile, which is almost 
certainly not maintained solely by high-Rossby number 
surface convection. Thus, even if there is a shallow re¬ 
turn flow above 0.977, this tells us little about the mean 
circulation at the base of the convection zone, which is 
most relevant for solar and stellar dynamo models (e.g. 
"Charbonneau |201 0). 


7.3. The Role of the Tachocline 
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Much as the upper boundary layer of the CZ might 
influence the meridional circulation profile of the solar 
convection zone, so too might the lower boundary layer, 
where the t rans ition to a sub-adiabatic radiative zone 
occurs. In §5.3| we attributed the spin-up of the poles 
in the high Rossby number regime to the poleward ad- 
vection of angular momentum by the extended single-cell 
meridional circulation. The circulation was attributed to 
gyroscopic pumping induced by an inward angular mo¬ 
mentum transport but positive feedback from an equa- 
torward entropy gradient amplified it, quickly establish¬ 
ing an anti-solar ft profile. Balance was only achieved 
when the shear became large enough for subgrid-scale 
diffusion to suppress further growth. 

In actual stars where (microscopic) diffusion is negli¬ 
gible, the balance may be achieved quite differently. For 
instance, thermal coupling to the tachocline can induce 
a baroclinicity throughout the CZ that helps warm the 
poles and thus suppresses the formation of a counter¬ 
clockwise (clockwise) circulation cell in the northern 
(southern) hemis phere ( |Rempel|2005 Miesch et al.| 2006 


simulations by Fan & Fang (2014). By suppressing the 


Brun et ah]|2011| . This baroclinicity is primarily due to 


the - downward circulation at high-latitudes which bur¬ 
rows down into the steep sub-adiabatic stratification of 
the overshoot region and radiative interior. This down¬ 
ward advection of relatively low entropy fluid from the 
CZ induces a strong poleward entropy gradient in the 
tachocline that establishes an approximate thermal wind 
balance with the differential rotation as expressed in eq. 
(15). Convective heat transport may then transmit this 
strong poleward entropy gradient vertically, back up into 
the CZ where it could mitigate the equatorward entropy 
gradient that hel ps drive the runaway convective mode 
discussed in §5.3| 

In addition to maintaining conical ft profiles through¬ 
out the CZ, this mechanism of thermal coupling may also 
suppress the strong CCW meridional circulation cell that 
helps to regulate the transition between solar and anti¬ 
solar rotation regimes. In short, the same mechanism 
that maintains the conical ft profile in the Sun might 
also help maintain the fast equator and slow pole that 
defines the solar-like profile (Aft > 0), shifting the anti¬ 
solar transition (Aft > 0) toward lower Rossby number. 
The conical ft profile of the Sun may thus signify that 
is i s clo se to the solar/anti-solar transition, as discussed 
in O We are currently investigating these issues with 
convection simulations that include a subadiabatic over¬ 
shoot region and radiative zone, to be described in future 
papers. 

7.4. Magnetism and Mean-Flow Regimes 

Another issue that warrants further research is the role 
of magnetism. Magnetism may suppress shear and thus 
assume the role of SGS diffusi on in the limit of low mi¬ 
croscopic diffusivity (e.g. Nelson et al. 2013). By sup¬ 


pressing shear, the Lorentz force (and in particular the 
Maxwell stress component) tends to partially offsetMhe 
zonal components of the Reynolds stress (Brun et al. 
2004| |Nelson et al.||2013|). This has the effect of reduc- 
mg the net axial torque T in equation (16), which in 


turn can reduce the amplitude of the meridional flow. 
This might account for why the presence of the magnetic 
field was found to shift the solar/anti-solar transition to¬ 
ward lower Rossby number in recent global convection 


meridional flow, magnetism may inhibit the polar spin- 
up associated with the anti-solar regime. 

The inhomogeneous spatial distribution of magnetic 
stresses can also suppress or alter meridional flow pat¬ 
terns. To illustrate the potential consequences of this, 
consider the rapidly-rotating regime, where the cylindri- 
cally outward angular momentum transport at low lat¬ 
itudes maintains a multi-cell meridional flow as shown 
in Fig. The region of divergence near the equator 
(red) is cylindrically aligned but its amplitude decreases 
with increasing distance from the equatorial plane. It is 
this axial variation in the Reynolds stress (parallel to the 
rotation axis) that helps establish the multi-cellular cir¬ 
culation pr ofile. In recent simulatio ns of convection in F- 
type stars, Augustson et al. (2013) find that the Lorentz 
force exerts a retrograde torque (flux divergence) at mid¬ 
latitudes that supplements the lower-latitude Reynolds 
stress, and diminishes its axial variation. The resulting 
meridional circulation at the equator in the lower CZ 
shifts from inward to outward, though the overall pro¬ 
file remains multi-celled. The presence of magnetism is 
thus likely to effect the detailed meridional flow profiles 
and the precise value of the transitional Rossby number 
between mean flow regimes, but it should not change 
our principle conclusion that rapid rotators should ex¬ 
hibit multi-cell meridional flow profiles while slow rota¬ 
tors should exhibit single-cell profiles. 

8. SUMMARY 

In this paper we report on a series of solar convec¬ 
tion simulations spanning a range of rotational influ¬ 
ences, as quan tified by the Rossby number Ro. As in 
previous work ([Gastine et al.||2013l |2014l Guerrero et al. 
2013a Kapyla et al. 2014), we find two distinct dynam- 
ical regimes of differential rotation (DR). One regime is 
characterized by solar-like DR (Aft > 0) at low Rossby 
number (Ro) and the other by anti-solar DR (Aft < 0) at 
high Ro, where A ft = ft eq — ft po i e is the pole-to-equator 
rotational contrast. Our transitional value of Ro ~ 0.16 
is lower than that reported in previous studies but we at¬ 
tribute this to a difference in how Ro is defined. The most 
novel aspect of the present work is a focus on the merid¬ 
ional circulation (MC) in the two dynamical regimes, in¬ 
cluding how it is established and maintained, its resulting 
cellular structure, and its role in defining the transition 
between regimes. 

We find a clear link between the DR and MC, with 
multi-cell MC profiles in the solar-like regime and single¬ 
cell MC profiles in the anti-solar regime. We attribute 
this to the convective angular momentum transport (i.e. 
the zonal component of the convective Reynolds stress, 
F^rs), which both establishes the DR contrast A ft and 
induces meridional flows through the mechanism of gyro¬ 
scopic pumping. In the solar-like regime, the positive A ft 
is maintained by a cylindrically outward transport (F#£ 
directed away from the rotation axis) at low latitudes. 
This induces a multi-cell MC profile outside the tangent 
cylinder, with two or more cells spanning the convection 
zone (CZ) at low latitudes, approximately aligned with 
the rotation axis. At high latitudes radially inward angu¬ 
lar momentum transport establishes a single circulation 
cell, with poleward flow in the upper CZ and equator- 
ward flow in the lower CZ. 
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As the rotational influence decreases, the qualitative 
nature of Frs changes. We demonstrate this change by 
introducing a transport potential x that highlights the 
orientation of Frs- In particular, the high-latitude re¬ 
gion of inward angular momentum transport extends to 
progressively lower latitudes as Ro is increased, while 
the low-latitude transport becomes equatorward. This 
establishes a single-celled MC profile throughout the CZ 
which is enhanced by baroclinic forcing. Poleward an¬ 
gular momentum transport by this induced circulation 
is mainly responsible for reversing the sign of AD and 
pushing the system into the anti-solar regime. 

Thus, though Frs ultimately regulates the transition 
between mean flow regimes, it is the angular momen¬ 
tum transport by the MC, and not the Reynolds stress 
itself, that leads to the spin-up of the polar regions in 
anti-solar states. This implies that additional factors im¬ 
pacting the efficiency of angular momentum transport by 
the MC, outside the purview of the simulations presented 
here, may play an important role in determining whether 
a particular star may attain an anti-solar state. These 
factors include magnetism and subtle boundary layer in¬ 
fluences involving the tachocline and near-surface shear 
layer. The dynamics near the transition are subtle. This 
is clear from the transitional case BO.5, which exhibits 
characteristics from both regimes, namely a solar-like DR 
and a single-celled MC. Though we did not investigate it 
explicitly here, we expect our simulations (and perhaps 
real stars) may also exhibit bistabi li ty ne ar the transi¬ 
tion as repo rted by |Gastine et al.| (|2014|) and |Kapyla 
et~aT| <f20l4|>. 


ent than in mean-field models where meridional flows 
are established by deviations from TWB du e to merid- 


We emphasize that Frs is the dominant factor that 
determines the MC profiles in our simulations through 
the mechanism of gyroscopic pumping. This mechanism 
can operate even if thermal wind balance (TWB), eq. 
(15) is satisfied everywhere. Thus, it is notably differ- 


ional Reynolds or Maxwell stresses (Kitchatinov 2012 
Dikpati 2014). These deviations are often localized in 


the boundary layers, which can have a disproportionate 
influence on the global MC profile. Real stars likely have 
significant contributions from both mechanisms, further 
increasing the challenge of predicting the MC profiles for 
individual stars. This is particularly the case for our own 
Sun, which may lie close to the transition between mean 
flow regimes. Indeed, the shallow return flow inferred 
from recent helioseismic inversions and photospheric fea¬ 
ture tracking may be a boundary-layer effect, reflect¬ 
ing the pen etrat ion depth of surface-driven convective 
plumes (sec. 7.2). 

Nevertheless, the results presented here suggest that at 
the extremes of the rotation spectrum, the convective an¬ 
gular momentum transport is likely the dominant driver 
of DR and meridional circulation. While the precise lo¬ 
cation of the transitional Ro may depend on a number of 
factors, solar-like stars rotating significantly faster than 
the Sun should thus have cylindrical, solar-like DR pro¬ 
files (AD > 0) and multi-celled MC profiles while those 
rotating significantly slower than the Sun should have 
anti-solar DR profiles (AD < 0) and single-celled MC 
profiles. 
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